Wstęp

Polski kodeks karny powstając w 1997 r. nie nadał jednoznacznie określonej definicji przestępstwa. Jednakże na podstawie przepisów możliwe jest jej wyodrębnienie: za przestępstwo rozumie się czyn pod groźbą kary, bezprawny, zawiniony i społecznie szkodliwy (Sobański, 2019). Badaczy, a dokładniej ekonomistów zastanawiało, co tak naprawdę determinuje bezprawne czyny. Według teorii, przestępczość i sprawiedliwość oparte jest na racjonalnych decyzjach jednostek, będące w ich interesie. Pierwowzorem ekonomicznej teorii przestępczości jest artykuł G. Beckera “Crime and Punishment: an economic approach”, który ukazał się w 1968 r. Przedstawia on rozumienie przestępstwa w ramach nurtu ekonomicznej analizy prawa. Jednostka dopuszcza się nielegalnego czynu, kiedy oczekiwana użyteczność z jego popełnienia jest większa od użyteczności, jaką mógłby osiągnąć poświęcając czas na inne aktywności. Podczas podejmowania decyzji o dokonaniu przestępstwa, osoba kieruje się prostym rachunkiem zysku i strat. Dopuści się go, jeżeli korzyści będą wyższe niż koszty (G.Becker, 1968). W związku z istniejącą teorią ekonomiczną przeprowadzono badanie determinantów przestępczości w powiatach. Opierając analizę na modelowaniu przestrzennym, a dokładniej modelu GWR (Geographically Weighted Regression)

Przegląd literatury

Analiza determinant przestępczości w literaturze naukowej wykazuje wieloaspektowe podejście badawcze, obejmujące czynniki ekonomiczne, społeczne, demograficzne oraz geograficzne. W podejściu ekonomicznym szczególnie istotne jest nawiązanie do klasycznej teorii G. Beckera, który wskazuje na bezrobocie jako czynnik redukujący koszt utraconych możliwości dla potencjalnych sprawców, a zarazem podkreśla znaczenie wysokości dochodów z legalnej działalności oraz ryzyka skazania dla kalkulacji korzyści i strat wynikających z działań nielegalnych (Becker, 1968). W tym kontekście korelacja pomiędzy przestępczością a rynkiem pracy pozostaje kluczowym obszarem badań.

Z kolei K. Kądziołka koncentruje się na analizie długookresowego bezrobocia, wykazując, że w Polsce problem ten jest bardziej dotkliwy dla mężczyzn, którzy są nad reprezentowani w strukturze sprawców przestępstw przeciwko mieniu. W jej badaniach szczególnie istotne okazały się wskaźnik urbanizacji oraz gęstość zaludnienia jako zmienne powiązane ze współczynnikiem przestępczości (Kądziołka, 2014). W podobnym tonie wypowiada się W. Florczak, wskazując na potrzebę precyzyjnego doboru regresorów w modelach przestępczości oraz uwzględnienia czynników odstraszających, takich jak wykrywalność przestępstw i surowość kar, które często są pomijane w analizach (Florczak, 2012).

W ujęciu społecznym istotnym elementem jest relacja pomiędzy dochodami a poziomem przestępczości. J.J. Sztaudynger i M. Sztaudynger zauważają, że wzrost przeciętnych dochodów może zarówno zmniejszać przestępczość poprzez zwiększenie kosztu utraconych korzyści, jak i ją intensyfikować poprzez zwiększenie atrakcyjności potencjalnych zysków z działań nielegalnych. Podkreślają oni również istotność spożycia alkoholu jako czynnika zaburzającego zdolność podejmowania racjonalnych decyzji i zwiększającego prawdopodobieństwo działań impulsywnych (Sztaudynger, Sztaudynger, 2003).

Perspektywa geograficzna badań nad przestępczością wskazuje na istotną rolę urbanizacji i koncentracji ludności. I. Bąk w swojej analizie podkreśla, że najwyższe wskaźniki przestępczości dotyczą dużych miast, natomiast analiza przestrzenna pozwala wyodrębnić różnice pomiędzy wschodnią a zachodnią częścią Polski. Wskazuje ona także na rozbieżności między oficjalnymi statystykami policyjnymi a percepcją zagrożenia przestępczością w społeczeństwie (Bąk, 2015).

Metody ekonometryczne stosowane w badaniach wskazują na złożoność interakcji między czynnikami ekonomicznymi a przestępczością. P. Bieniek, S. Cichocki i M. Szczepaniec przeprowadzili analizę wskazującą na dodatnią korelację pomiędzy przestępczością a wielkością miasta oraz bliskością do zachodniej granicy Polski, przy czym czynniki takie jak edukacja, wydatki na bezpieczeństwo i zużycie energii elektrycznej wykazywały odwrotną korelację z liczbą przestępstw (Bieniek, Cichocki, Szczepaniec, 2012). Podobne rezultaty uzyskano w badaniach S. Mordwy dotyczących przestrzennej dystrybucji kradzieży, gdzie wykazano, że czynniki demograficzne takie jak gęstość zaludnienia oraz struktura płciowa wpływają na intensywność przestępstw przeciwko mieniu (Mordwa, 2011). Natomiast metody przestrzenne, jak GWPR podkreślają złożone relację między nierównościami dochodów a przestępczością przeciwko mieniu, a sprzeczne wyniki dotyczące relacji między stopami bezrobocia, czy czynnikami niskich dochodów a wskaźnikami włamań. Ponadto w badanach przejawia się złożona relacja między imigracją a przestępczością przeciwko mieniu (Traves, 2021).

Dodatkowo, K. Kądziołka przeprowadziła analizę wpływu czynników ekonomicznych, społecznych i demograficznych na przestępczość w Polsce, wykorzystując metodę drzew regresji. Badanie wykazało, że powiaty zachodniej części kraju charakteryzują się ponadprzeciętnym poziomem przestępczości, a kluczowe determinanty to wskaźnik urbanizacji, gęstość zaludnienia, odsetek gospodarstw jednoosobowych oraz natężenie przestępstw w sąsiednich powiatach. Zastosowanie metody lasu losowego pozwoliło na redukcję skokowości danych i umożliwiło przeprowadzenie analizy na poziomie globalnym. Wyniki te dostarczają istotnych wskazówek dla polityki prewencyjnej oraz optymalizacji strategii bezpieczeństwa (Kądziołka, 2016).

W kontekście analizy systemowej przestępczości G. Dubiec dokonał zastosowania metod wielowymiarowej analizy porównawczej, wskazując na podział determinant na czynniki ekonomiczne, społeczne, demograficzne i geograficzne. Kluczowe zmienne obejmowały m.in. poziom dochodów, wydatki socjalne, liczbę eksmisji oraz wskaźnik urbanizacji. Zastosowane techniki statystyczne pozwoliły na wskazanie regionów o podwyższonej przestępczości oraz wypracowanie rekomendacji dla polityki prewencyjnej (Dubiec, 2020). Podobne wnioski wynikają z badań Ł. Gadomskiego i M. Malczewskiego, którzy, opierając się na modelach panelowych, wykazali silny związek pomiędzy warunkami życiowymi a poziomem przestępczości. Ich wyniki wskazują, że wzrost poziomu życia o 1% powoduje spadek przestępczości o 0,45% (Gadomski, Malczewski, 2020).

Metoda Badawcza

Jak zauważa Tavares (2021) badania nad wpływem środowiska na przestępczość coraz częściej wykorzystują narzędzia geoinformacyjne, eksplorację danych oraz wielowymiarową analizę statystyczną. Efekty przestrzenne, takie jak autokorelacja i heterogeniczność przestrzenna, wymagają zastosowania modeli uwzględniających te specyficzne zależności. W kontekście analizy przestępczości szeroko stosowane są modele opóźnienia przestrzennego i błędu przestrzennego, a także podejście autoregresji przestrzennej Cliffa-Orda. Coraz większe znaczenie zyskuje również regresja ważona geograficznie (GWR), która umożliwia uwzględnienie lokalnych zmienności w analizowanych zależnościach. W badaniach nad przestępczością stosuje się także modele regresji Poissona oraz regresji ujemnej dwumianowej, które pozwalają na modelowanie danych o charakterze zliczeniowym. Geograficzne odmiany tych metod, takie jak GWPR i GWNBR, stanowią cenne narzędzie do analizy przestrzennych wzorców przestępczości i ich determinant (Tavares, 2021). Tym samym, analizę w tej pracy również oparto na metodzie GWR.

Na podstawie literatury sformułowano następujące hipotezy:

H1: Zgodnie z ekonomiczną teorią przestępczości; wyższa wartość dochodu powinna obniżać przestępczości, ponieważ zwiększają się koszty utraconych możliwości. W tym przypadku założono, że zarobki zmniejszają przestępczość - wyższe średnie wynagrodzenie w powiatach będzie obniżać ilość popełnionych przestępstw. (K. Kądziołka, Wpływ wybranych czynników o charakterze społeczno-ekonomicznym na przestępczość przeciwko mieniu w Polsce, 2014, s. 22).

H2: Silny związek przestępczości z bezrobociem. Brak pracy, a w szczególności długotrwałe bezrobocie budzi poczucie wykluczenia, czy też braku nadziei na uzyskanie dochodowej pracy. Zwiększa to motywację do działań nielegalnych (J.J. Sztaudynger, M. Sztaudynger, Ekonometryczne modele przestępczości, 2003).

H3: Występowanie efektu odstraszania ogólnego, który definiuje zmienna wskaźnika wykrywalności przestępstw. Skuteczność wymiaru sprawiedliwości będzie zniechęcała podmioty do popełniania przestępstw.

H4: Wysoki współczynnik rozwodów i separacji może mieć negatywny wpływ na ogólny dobrobyt społeczeństwa. Może to zwiększyć prawdopodobieństwo wyższej przestępczości na danym terenie (P. Bieniek, S. Cichocki, M. Szczepaniec, Wyniki ekonomiczne a poziom przestępczości - badanie ekonometryczne, 2012).

H5: Wskaźnik urbanizacji jest pozytywnym czynnikiem podaży przestępczości. Ośrodki miejskie charakteryzują się większą ilością odchyleń społecznych dotyczących przestrzegania prawa. Spowodowane jest to znaczącym ułatwieniem popełnienia przestępstwa, miasta odznaczają się większymi możliwościami w osiągnięciu zysku.

H6: Kazus gęstości zaludnienia jest podobny, co do wskaźnika urbanizacji; znacznie zwiększa możliwość dokonania przestępstwa. Regiony o większej gęstości zaludnienia charakteryzują się większą stopą zwrotu niż wsie czy małe ośrodki miejskie.

H7: Wyższa wartości dochodu rozporządzalnego regionu obniża ilość przestępstw. Im wyższe odchylenie od progu ubóstwa tym obywatele o wyższych dochodach otrzymują mniej bodźców na rzecz łamania prawa. Jeśli dochody faktycznie należą do wysokich to jest to zniechęcenie do popełniania przestępstw ze względu na rosnące koszty utraconych możliwości.

H8: Bardziej wykształcone społeczeństwo powinno się charakteryzować niższym wskaźnikiem przestępczości. Lepsza edukacja może zwiększać świadomość prawnospołeczną i zwiększać atrakcyjność na rynku pracy. Bardzo ważne w tym przypadku jest jak wysoki poziom jest osiągany społecznie.  Niektóre źródła zwracają uwagę, że przestępczość wśród osób gorzej wykształconych jest związana z ich wykorzystywaniem w nielegalnych działalnościach - często są wykorzystywane jako „słupy” w przestępstwach polegających na praniu brudnych pieniędzy (Ł. Gadomski, M. Malaczewski, Zastosowanie dynamicznych modeli panelowych w analizie determinant przestępczości w Polsce w latach 2005–2016, 2020). Analogicznie, im więcej osób z niskim wykształceniem, tym wyższa powinna być przestępczość.

H9: Kolejnym czynnikiem mający wpływ na całkowitą wartość przestępczości jest współczynnik feminizacji. Wyższy odsetek kobiet w regionie wpływa na obniżenie ilości przestępstw, ponieważ kobiety charakteryzują się niższą skłonnością do popełnienia przestępstwa.

H10: Co do edukacji kluczowe również są wydatki na edukacje na jednego obywatela. Im wyższe wartości kwotowe są osiągane tym możemy założyć, iż szanse na popełnienie przestępstwa będą malały. Informuję nas to o tym, iż społeczeństwo dalej się kształci w jakimś stopniu.

H11: W przypadku imigracji, czynnik ten ma za zadanie uwzględnić kontekst społeczno-ekonomiczny, a głównie środowisko występowania przestępstw. Powiaty z większa liczbą imigrantów, będące bardziej zróżnicowane etniczne będą wykazywały zwiększoną liczbę nielegalnych czynów (Wang et al. 2019).

H12: Jednoosobowe gospodarstwa ze względu na mniejszy poziom ochrony niż gospodarstwa wieloosobowe są bardziej narażone na włamania, czy też kradzieże. Tym samym, wyższy odsetek takich gospodarstw będzie zwiększał dostępność potencjalnych ofiar, tym samym przyczyniał się do zwiększonej liczby przestępstw.

Dane

W analizie wykorzystano dane pochodzące z Banku Danych Lokalnych, a dokładniej Narodowego Spisu Powszechnego, który miał miejsce w 2021 r. Uzyskując informację dla gospodarstw domowych ze wszystkich powiatów w Polsce, których liczba wynosi 380. Badana zmienna to logarytm ilości przestępstw na 10 000 mieszkańców. Zmienne objaśniające, to zarówno mierniki ekonomiczne, demograficzne oraz geograficzne, pełna lista to:

Zmienna Opis
log_przestepstw Logarytm ilości przestępstw na 10 000 mieszkańców popełnionych na danym obszarze terytorialnym.
DochodyNaMieszkanca Dochód budżetów powiatów przypadający na 1 mieszkańca.
Feminizacja Stosunek liczby kobiet przypadających na 100 mężczyzn.
Imigranci Liczba imigrantów przebywających w danym regionie.
Wykrywalnosc Stosunek procentowy liczby przestępstw wykrytych do ogólnej liczby przestępstw stwierdzonych.
rozwody_na10000 Średnia liczba rozwodów przypadająca na 10 000 mieszkańców na danym terenie.
bezrobocie Stosunek liczby osób zarejestrowanych jako bezrobotne do liczby ludności aktywnej zawodowo.
wydatki_bezp Wysokość wydatków na bezpieczeństwo w regionach.
wydatki_oswiata Wysokość wydatków na oświatę i edukację w regionach.
wydatki_pomocspol Wysokość wydatków na pomoc społeczną w regionach.
wynagrodzenie_relacja Wysokość średnich zarobków w regionach w relacji do średniej krajowej =100.
urban Procentowy udział mieszkańców miast w ogólnej liczbie ludności.
zalud Gęstość zaludnienia, miara przyporządkowująca populację (ilość osób) do danej powierzchni (km2).
wyk1 Procent osób z wykształceniem poniżej średniego.
gosp_jedno Odsetek jednoosbowych gospodarstw domowych.

Wizualizacja danych

par(mfrow = c(1, 2))
zmienna<-DaneFinal$przestepstwa_na10000
maxy<-95
breaks<-c(5, 15, 25, 35, 45, 55, 65, 75, 85, 95)  

nclr<-8
plotclr<-brewer.pal(nclr, "Reds") # {RColorBrewer} 
fillRed<-colorRampPalette(plotclr) # {grDevices}  
colcode<-fillRed(maxy)[round(zmienna) + 1] #fillRed 

par(mar=c(4,4,4,4))
plot(st_geometry(POW), col = colcode, lty = 0, border = "gray")  
plot(st_geometry(WOJ), add=TRUE, lwd=1, border = "gray60")  
colorlegend(posy=c(0.05,0.9), posx=c(0.1,0.12), col=fillRed(maxy), zlim=c(0, maxy), zval=breaks, main.cex=0.9)

Mapa powyżej przedstawia rozkład przestępstw w Polsce. Średnia liczba przestępstw w powiatach w 2021 r. wynosiła 20 na 10 000 mieszkańców, jednak w poszczególnych jednostkach terytorialnych wartości te znacząco odbiegały od średniej krajowej. Najwyższym wskaźnikiem przestępczości odznaczał się Nowy Sącz, gdzie na 10 000 mieszkańców przypadało 91 przestępstw, co stanowi najwyższą wartość w skali kraju. Kolejne pozycje zajmowały powiaty sieradzki (91), sierpecki (71), m. Katowice (69) oraz m. Krosno (57). Dominacja dużych miast oraz obszarów o wysokim stopniu zurbanizowania może wynikać z większej liczby potencjalnych sprawców i ofiar. Z kolei najniższe wartości wskaźnika przestępczości odnotowano w powiatach strzyżowskim (8), rzeszowskim (8), proszowickim (7), przemyskim (7) oraz grudziądzkim (7). Są to w większości regiony o niskiej gęstości zaludnienia oraz przewadze terenów wiejskich, co może sprzyjać niższej skali czynów zabronionych, wynikających z ograniczonej liczby potencjalnych ofiar przestępstw.

library(classInt)
library(RColorBrewer)

par(mfrow = c(2, 2), mar = c(3, 3, 3, 6)) 

vars <- list(
  "DochodyNaMieszkanca" = c(500, 1500, 4500, 9500, 12500, 15500, 20000),
  "wydatki_bezp" = c(10000, 5000000, 10000000, 50000000, 100000000, 200000000, 330000000),
  "wydatki_oswiata" = c(0, 200, 400, 600, 1000, 2000, 3500),
  "wydatki_pomocspol" = c(300000, 5000000, 10000000, 50000000, 100000000, 500000000, 650000000)
)

for (var in names(vars)) {
  zmienna <- DaneFinal[[var]]
  breaks <- vars[[var]]
  przedzialy <- length(breaks) - 1  
  
  kolory <- brewer.pal(przedzialy, "BuPu")  
  klasy <- classIntervals(zmienna, przedzialy, style="fixed", fixedBreaks=breaks)  
  tabela.kolorow <- findColours(klasy, kolory)  
  
  plot(st_geometry(POW), col=tabela.kolorow, lty=0, border="gray")  
  plot(st_geometry(WOJ), add=TRUE, lwd=1, border="black")  
  
  
  legend("topright", legend=names(attr(tabela.kolorow, "table")), 
         fill=attr(tabela.kolorow, "palette"), cex=0.8, bty="n", 
         xpd=TRUE, inset=c(-0.2, 0))  
  

  title(main=var, cex.main=1.2, line=1)  
}

Dochody budżetów powiatów przypadające na jednego mieszkańca wykazują silne zróżnicowanie. Średnia wartość dla wszystkich powiatów wynosi 2710,8 zł, a najwyższe dochody odnotowano w miastach na prawach powiatu, w tym w m. Świnoujście (19 462 zł), m. Sopot (13 020 zł) i m. Warszawa (11 608 zł). Najniższe wartości osiągnęły natomiast powiaty o niższym poziomie turystycznym i uprzemysłowienia, jak rybnicki (841 zł), skierniewicki (846 zł) i opolski (857 zł).

W przypadku wydatków na bezpieczeństwo publiczne, oświatę czy pomoc społeczną można oczekiwać, że powiaty z wyższymi budżetami będą rozdysponowywać więcej pieniędzy pomiędzy mieszkańcami. Średni poziom wydatków na bezpieczeństwo w skali kraju wynosi 10,6 mln zł, jednak pomiędzy regionami występują znaczące dysproporcje. Najwyższe nakłady odnotowywane są w największych ośrodkach miejskich, jak Warszawa (329 mln zł), Kraków (150 mln zł) i Łódź (116 mln zł). Z kolei powiaty o najmniejszych wydatkach, to skierniewicki (12 tys. zł), łomżyński (23 tys. zł) i krośnieński (24 tys. zł).

Podobnie zróżnicowane są nakłady na oświatę i edukację, gdzie średnia wartość na mieszkańca wynosi 794,4 zł. Jednakże w niektórych powiatach te wartości są kilkukrotnie wyższe. Najwyższe nakłady ponoszone są w m. Siedlcach (3 399 zł), m. Zamość (3 269 zł) i m. Ostrołęka (3 194 zł). W przeciwieństwie do tego, powiaty łomżyński (0,8 zł), koszaliński (43 zł) i zamojski (61 zł) wydają na edukację nawet kilkudziesięciokrotnie mniej od średniej krajowej. Można przypuszczać, że w wyniku tak niskich nakładów na rozwój oświaty w powiatach tych będzie występował wyższy odsetek osób z wykształceniem poniżej średniego.

Wydatki na pomoc społeczną wykazują najwyższe zróżnicowanie. Średnia wartość wynosi 23,2 mln zł. Największe budżety na pomoc społeczną jak w przypadku reszty wydatków posiadają największe aglomerację, w tym m. Warszwa (646 mln zł), m. Kraków (414 mln zł) i m. Łódź (298 mln zł). Znacznie niższe wydatki występują w powiecie skierniewickim (377 tys. zł), łosickim (405 tys. zł) i grudziądzkim (444 tys. zł).

library(classInt)
library(RColorBrewer)
par(mfrow = c(1, 3), mar = c(3, 3, 2, 3))  # Układ 1x3, dostosowane marginesy

vars <- list(
  "Feminizacja" = c(95, 100, 105, 110, 115, 120),
  "Imigranci" = c(0, 500, 1000, 5000, 10000, 50000, 170000),
  "rozwody_na10000" = c(4, 10, 13, 15, 18, 20, 28)
)

for (var in names(vars)) {
  zmienna <- DaneFinal[[var]]
  breaks <- vars[[var]]
  przedzialy <- length(breaks) - 1  

  kolory <- brewer.pal(przedzialy, "BuPu")  
  klasy <- classIntervals(zmienna, przedzialy, style="fixed", fixedBreaks=breaks)  
  tabela.kolorow <- findColours(klasy, kolory)  

  plot(st_geometry(POW), col=tabela.kolorow, lty=0, border="gray")  
  plot(st_geometry(WOJ), add=TRUE, lwd=1, border="black")  

  legend("bottomleft", legend=names(attr(tabela.kolorow, "table")), 
         fill=attr(tabela.kolorow, "palette"), cex=1, bty="n")

  title(main=var, cex.main=1.0, line=-1)  
}

Zmienne jak wskaźnik feminizacji, liczba imigrantów i ilość rozwodów odzwierciedlają kontekst społeczny w danych powiatach. Feminizacja określa stosunek liczby kobiet do liczby mężczyzn w danym regionie, średnia dla kraju wynosi 105. Najwyższe wartości odnotowano w dużych miastach, jak m. Łódź (119), m. Sopot (119) oraz m. Warszawa (117). Wysoki wskaźnik w miastach może wynikać z większej migracji kobiet do tych regionów w celu podjęcia pracy, czy edukacji. Z kolei najniższe wskaźniki występują w powiatach wiejskich, takich jak suwalski (95), ostrołęcki (97) i łomżyński (98).

Średnia liczba imigrantów przebywających w powiatach wynosi 3 733, jednak między regionami występują duże różnice. Największe skupiska imigrantów znajdują się w głównych ośrodkach miejskich, zwłasza w Warszawie (166 691). Następnie wysoka liczba imigrantów odnotowywana jest w Wrocławiu (69 099) i Krakowie (54 630). Są to miasta wykazujące się atrakcyjnością ekonomiczną, a zwłaszcza edukacyjną. Z kolei powiaty o najniższej liczbie imigrantów, to węgorzewski (78), kolbuszowski (115) i kolneński (128).

Na 10 000 mieszkańców średnio występuje 15 rozwodów, jednakże w wybranych powiatach wskaźnik jest wyższy. Tym razem nie są to główne aglomeracje, a powiaty takie jak m. Konin (27,3), m. Włocławek (24,5) i m. Biała Podlaska (23,5), jednakże nadal są to ośrodki miejskie. Najniższe wskaźniki odnotowano w powiatach proszowickim (4,5), suskim (4,6) i wadowickim (6,2), które odznaczają się niskimi wskaźnikami urbanizacji. Są to regiony wiejskie, stąd niski wskaźnik rozwodów może wynikać z silniejszego wpływu tradycyjnych wartości na życie rodzinne, bądź niskiego poziomu dochodów.

library(classInt)
library(RColorBrewer)
par(mfrow = c(1, 3), mar = c(3, 3, 2, 3))  

vars <- list(
  "urban" = c(0, 0.2, 0.3, 0.4, 0.5, 0.7, 1),
  "zalud" = c(10, 50, 100, 200, 500, 1000, 3700),
  "gosp_jedno" = c(0.1, 0.15, 0.17, 0.19, 0.21, 0.25, 0.36)
)

for (var in names(vars)) {
  zmienna <- DaneFinal[[var]]
  breaks <- vars[[var]]
  przedzialy <- length(breaks) - 1  

  kolory <- brewer.pal(przedzialy, "BuPu")  
  klasy <- classIntervals(zmienna, przedzialy, style="fixed", fixedBreaks=breaks)  
  tabela.kolorow <- findColours(klasy, kolory)  

  plot(st_geometry(POW), col=tabela.kolorow, lty=0, border="gray")  
  plot(st_geometry(WOJ), add=TRUE, lwd=1, border="black")  

  legend("bottomleft", legend=names(attr(tabela.kolorow, "table")), 
         fill=attr(tabela.kolorow, "palette"), cex=1, bty="n")

  title(main=var, cex.main=1.0, line=-1)  
}

Średni procentowy udział ludności miejskiej w powiatach wynosi 50,4%. Najwyższy wskaźnik urbanizacji (pomijając miasta na prawach powiatu, gdzie wskaźnik ten wynosi 1) odnotowano w powiatach mikołowskim (85%), dzierżoniowskim (78%) i bieruńsko-lędzińskim (75%). Są to też powiaty położone w województwie śląskim, w bliskim sąsiedztwie do dużych aglomeracji jak Katowice czy Wrocław. Z kolei najniższy poziom urbanizacji występuje w powiatach skierniewickim (0%), suwalskim (0%) i przemyskim (1,1%). Są to obszary o silnie wiejskim charakterze, położone z daleka od dużych ośrodków miejskich, czy blisko wschodniej granicy Polski.

Gęstość zaludnienia w powiatach wynosi średnio 350 osób/km2, najwyższe wskaźniki gęstości odnotowano w największym polskim mieście Warszawa (3 606 osób/km2). Następnie są to Świętochłowice (3 501) i Chorzów (3 106), położone w województwie śląskim, który charakteryzuje się wysoką urbanizacją i gęstością zaludnienia. Z kolei najniższe gęstości odnotowano w powiecie bieszczadzkim (17 osób/km2), sejeńskim (21) i hajnowskim (24).

Odsetek jednoosobowych gospodarstw domowych w powiatach wynosi średnio 19,4%. Najwyższy udział widoczny jest w Warszawie (35%), Sopocie (34%) i Łodzi (34%), czyli ponownie w dużych aglomeracjach, gdzie struktura społeczna miast spzryja modelowi życia indywidualnego. Z kolei najniższe odsetki jednoosboowych gospodarstw odnotowano w powiatach leszczyńskim (10,6%), nowosądeckim (10,7%) i konińskim (10,8%).

library(classInt)
library(RColorBrewer)
par(mfrow = c(1, 3), mar = c(3, 3, 2, 3))  # Układ 1x3, dostosowane marginesy

vars <- list(
  "wynagrodzenie_relacja" = c(70, 80, 90, 100, 110, 120, 170),
  "wyk1" = c(0.1, 0.3, 0.35, 0.4, 0.42, 0.45, 0.49),
  "bezrobocie" = c(1, 3, 4, 5, 6, 8, 12)
)

for (var in names(vars)) {
  zmienna <- DaneFinal[[var]]
  breaks <- vars[[var]]
  przedzialy <- length(breaks) - 1  

  kolory <- brewer.pal(przedzialy, "BuPu")  
  klasy <- classIntervals(zmienna, przedzialy, style="fixed", fixedBreaks=breaks)  
  tabela.kolorow <- findColours(klasy, kolory)  

  plot(st_geometry(POW), col=tabela.kolorow, lty=0, border="gray")  
  plot(st_geometry(WOJ), add=TRUE, lwd=1, border="black")  

  legend("bottomleft", legend=names(attr(tabela.kolorow, "table")), 
         fill=attr(tabela.kolorow, "palette"), cex=1, bty="n")

  title(main=var, cex.main=1.0, line=-1)  
}

Wysokość wynagrodzeń w powiatach średnio wynosi 86,65% średniej krajowej, co oznacza, że przeciętne wynagrodzenia poza największymi ośrodkami miejskimi są niższe od średniej krajowej. Najwyższe wynagrodzenia odnotowano w powiecie lublińskim (167,9), m. Jastrzębie-Zdrój (149,8) i m. Warszawa (128,1). Wysokie wynagrodzenia w tych regionach wynikają z obecności dużych przedsiębiorstw, jak KGM w Lublinie, czy przypadku Jastrzębia-Zdrój rozwiniętego górnictwa węgla kamiennego. Najniższe wynagrodzenia występują natomiast w powiecie kępińskim (70,7), krotoszyńskim (71,4) i złotowskim (71,7), w których brak jest dużych zakładów pracy.

Średnia wartość udziału osób z wykształceniem poniżej średniego w powiatach wynosi 39,4%, co oznacza, że w wielu regionach ponad 1/3 mieszkańców nie ukończyła szkoły średniej. Najwyższy odsetek tych osób występuje w powiecie nowomiejskim (48,8%), sępoleńskim (48,5%) i rypińskim (48,3%). Z kolei najniższy udział osób z wykształceniem poniżej średniego odnotowano w dużych miastach, jak Warszawa (17,6%), Rzeszów (20,5%) i Kraków (21,6%). Miasta te oferują szersze możliwości edukacyjne oraz obecność uczelni wyższych, co sprzyja rozwojowi edukacji.

Średnia stopa bezrobocia w Polsce wynosi 4,58%, jednak wartości te różnią się między regionami. Najwyższe bezrobocie występuje w powiatach szydłowskim (11,8%), m. Przemyśl (10,3%) i kętrzyńskim (10,3%). Najniższe stopy odnotowano w powiecie kepińskim (1,5%), wolsztyńskim (1,8%) i grójeckim (1,9%).

Model OLS

model1 <- lm(log_przestepstwa ~ DochodyNaMieszkanca + Feminizacja + Imigranci + wykrywalnosc +
              rozwody_na10000 + bezrobocie + wydatki_bezp + wydatki_oswiata +
              wydatki_pomocspol + urban + zalud + wyk1 +
              wynagrodzenie_relacja + gosp_jedno, data = DaneFinal)

summary(model1)
## 
## Call:
## lm(formula = log_przestepstwa ~ DochodyNaMieszkanca + Feminizacja + 
##     Imigranci + wykrywalnosc + rozwody_na10000 + bezrobocie + 
##     wydatki_bezp + wydatki_oswiata + wydatki_pomocspol + urban + 
##     zalud + wyk1 + wynagrodzenie_relacja + gosp_jedno, data = DaneFinal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.67750 -0.17281 -0.02118  0.12558  1.43957 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            1.301e+00  9.901e-01   1.314  0.18961    
## DochodyNaMieszkanca    1.573e-06  1.667e-05   0.094  0.92488    
## Feminizacja           -2.081e-03  8.637e-03  -0.241  0.80971    
## Imigranci              6.208e-06  4.383e-06   1.416  0.15751    
## wykrywalnosc           1.171e-02  2.086e-03   5.616 3.88e-08 ***
## rozwody_na10000        1.532e-02  5.479e-03   2.796  0.00545 ** 
## bezrobocie            -3.328e-02  1.092e-02  -3.049  0.00247 ** 
## wydatki_bezp          -1.284e-09  3.341e-09  -0.384  0.70104    
## wydatki_oswiata       -5.427e-06  6.825e-05  -0.080  0.93667    
## wydatki_pomocspol     -9.844e-10  1.162e-09  -0.847  0.39745    
## urban                  4.718e-01  1.477e-01   3.194  0.00152 ** 
## zalud                 -4.022e-05  5.706e-05  -0.705  0.48136    
## wyk1                  -2.456e-01  4.368e-01  -0.562  0.57418    
## wynagrodzenie_relacja  1.544e-03  1.781e-03   0.867  0.38660    
## gosp_jedno             3.222e+00  7.236e-01   4.452 1.13e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2828 on 365 degrees of freedom
## Multiple R-squared:  0.4909, Adjusted R-squared:  0.4714 
## F-statistic: 25.14 on 14 and 365 DF,  p-value: < 2.2e-16
(vif_values <- car::vif(model1))
##   DochodyNaMieszkanca           Feminizacja             Imigranci 
##             10.151959              5.702086             10.277211 
##          wykrywalnosc       rozwody_na10000            bezrobocie 
##              1.400228              1.887690              1.393054 
##          wydatki_bezp       wydatki_oswiata     wydatki_pomocspol 
##             24.985491             14.968544             15.801309 
##                 urban                 zalud                  wyk1 
##              7.790650              6.002075              3.174370 
## wynagrodzenie_relacja            gosp_jedno 
##              1.707306              5.734176
(AIC_value <- AIC(model1))
## [1] 135.3075
bptest(model1) 
## 
##  studentized Breusch-Pagan test
## 
## data:  model1
## BP = 40.011, df = 14, p-value = 0.0002542

Początkowy model zawiera 14 zmiennych, gdzie tylko 5 z nich jest istotne statystycznie na poziomie istotności 0,05, są to: wskaźnik wykrywalności, rozwody, stopa bezrobocia, urbanizacja oraz odsetek gospodarstw jednodomowych. Wysokim wskaźnikiem VIF odznaczają się zmienne dot. dochodów na mieszkańca, wydatków na pomoc społeczną, wydatków na oświatę, wydatków na bezpieczeństwo oraz liczby imigrantów. Każda z tych zmiennych posiada współczynnik VIF powyżej 10, a nawet powyżej 20, co tym samym może zaburzać oszacowania w wyniku wysokiej współliniowości.

Przeprowadzając test Breuscha-Pagana na heteroskedastycznośc, uzyskano p-value równe 0, co prowadzi do odrzucenia hipotezy zerowej mówiącej o homoskedastyczności. Zbadano również statystykę AIC, która wyniosła 135.31, a następnie przystąpiono do próby poprawy modelu.

# Model
eq = log_przestepstwa~Feminizacja+log_Imigranci+wykrywalnosc+log_rozwody_na10000+
  bezrobocie+wynagrodzenie_relacja+log_wydatki_bezp+urban+log_zalud+log_gosp_jedno

model.lm=lm(eq, data=DaneFinal)
summary(model.lm)
## 
## Call:
## lm(formula = eq, data = DaneFinal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.57370 -0.17410 -0.03533  0.13146  1.37234 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            2.732619   1.050479   2.601  0.00966 ** 
## Feminizacja           -0.015005   0.008943  -1.678  0.09421 .  
## log_Imigranci          0.056976   0.018743   3.040  0.00254 ** 
## wykrywalnosc           0.012956   0.002025   6.397 4.80e-10 ***
## log_rozwody_na10000    0.234384   0.074178   3.160  0.00171 ** 
## bezrobocie            -0.013966   0.012435  -1.123  0.26211    
## wynagrodzenie_relacja  0.001334   0.001658   0.805  0.42154    
## log_wydatki_bezp       0.022630   0.013954   1.622  0.10572    
## urban                  0.272084   0.146135   1.862  0.06342 .  
## log_zalud              0.040681   0.026189   1.553  0.12119    
## log_gosp_jedno         0.601437   0.135913   4.425 1.27e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2762 on 369 degrees of freedom
## Multiple R-squared:  0.5091, Adjusted R-squared:  0.4958 
## F-statistic: 38.27 on 10 and 369 DF,  p-value: < 2.2e-16
(vif_values <- car::vif(model.lm))
##           Feminizacja         log_Imigranci          wykrywalnosc 
##              6.410258              2.429398              1.384499 
##   log_rozwody_na10000            bezrobocie wynagrodzenie_relacja 
##              1.888984              1.895398              1.551411 
##      log_wydatki_bezp                 urban             log_zalud 
##              1.875877              7.995582              5.284477 
##        log_gosp_jedno 
##              5.280520
(AIC_value <- AIC(model.lm))
## [1] 113.4606
bptest(model.lm) 
## 
##  studentized Breusch-Pagan test
## 
## data:  model.lm
## BP = 32.436, df = 10, p-value = 0.0003387

Uproszczony model zawiera 10 zmiennych, gdzie 4 z nich jest istotna statystycznie na poziomie istotności 0,05, a przyjmując poziom istotności 0,1 liczba istotnych zmiennych wzrasta do 6. Statystyka AIC poprawiła się znacząco, wynosi 113.46, dlatego można uznać, że model uproszczony jest lepszy od początkowego modelu. Wprowadzono logarytmy zmiennych, aby zbliżyć ich rozkład do normalnego i ułatwić interpretację ekonometryczną. Wyróżniającymi się zmiennymi podczas badania współczynnika VIF są: feminizacja, urbanizacja, logarytm zaludnienia i logarytm odsetka gospodarstw domowych, gdzie każda z tych zmiennych przyjmuje wartość VIF powyżej 5, co wskazuje na wysoką współliniowość. Jednakże ze względu na literaturę przedmiotu i ich wysoką istotność postanowiono pozostawić te zmienne w modelu. Test Breuscha Pagana (p-value = 0,00) ponownie wykazuje na występowanie heteroskedastyczności, co może sugerować potrzebę zastosowania innych metod ekonometrycznych, jak modelowanie przestrzenne.

Model OLS wskazuje ogólną zależność relacji, można zauważyć, że najbardziej odznaczającymi się zmiennymi są: logarytm rozwodów, logarytm odsetka gospodarstw jednoosbowoych, wykrwyalność oraz logarytm liczby imigrantów. Każda z tych zmiennych oddziałuje pozytywnie na liczbę stwierdzanych przestępstw, prowadząc do tego, że ich wzrost będzie przyczyniał się do wzrostu ilość nielegalnych czynów. Na granicy istotności znajdują się natomiast zmienne, jak: urbanizacja i feminizacja, których istotna zależność wykazywana jest w literaturze przedmiotu. Reszta zmiennych mimo nieistotności nie jest usuwana z modelu ponieważ mogą być przestrzenne heterogenicznie, dlatego potrzeba jest odwołania się do innych metod niż klasyczny model OLS.

plot(model.lm$fitted.values, residuals(model.lm), 
     xlab = "Wartości dopasowane", ylab = "Reszty", 
     main = "Wykres reszt")
abline(h = 0, col = "red")

Rozkład reszt modelu liniowego nie wykazuje jednoznacznego trendu ani kształtu, co sugeruje, że relacja między zmiennymi niezależnymi a zależna jest w miarę poprawnie dopasowana. Punkty są rozłozone mniej więcej równomiernie wokół osi poziomej, brak wyraźnego wzorca w wariancji reszt. Widać występowanie kilku wartości odstających. Jednakże warto przyjrzeć się przestrzennemu rozkładowi reszt.

res<-model.lm$residuals
brks<-c(min(res), mean(res)-sd(res), mean(res), mean(res)+sd(res), max(res))
cols<-c("steelblue4","lightskyblue","thistle1","plum3")
plot(pow, col=cols[findInterval(res,brks)])
title(main="Reszty w modelu MNK")
legend("bottomleft", legend=c("<mean-sd", "(mean-sd, mean)", "(mean, mean+sd)", ">mean+sd"), leglabs(brks), fill=cols, bty="n")

Analiza przestrzennego rozkładu reszt modelu umożliwia identyfikację regionalnych wzorców różnic między wartościami rzeczywistymi a przewidywaniami modelu. Obserwowane odchylenia wskazują na potencjalne niedoszacowania lub przeszacowania, co może wynikać z nieuwzględnionych czynników strukturalnych w modelu. Powiaty oznaczone kolorem ciemnoróżowym wykazują wartości rzeczywiste znacznie wyższe od przewidywanych, tam gdzie model przeszacował wartości widnieje kolor ciemnoniebieski. Wysokie reszty (ciemnoróżowe powiaty) dominują w południowo-wschodniej Polsce oraz części Mazowsza, a niskie reszty (ciemnoniebieskie obszary) we wschodniej i północnej części kraju. Możliwe jest istnienie silnych zależności przestrzennych, jednak aby je okreslić i poprawić jakość modelu należy przeprowadzić test Morana.

lm.morantest(model.lm, cont.listw)
## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = eq, data = DaneFinal)
## weights: cont.listw
## 
## Moran I statistic standard deviate = 3.1288, p-value = 0.0008777
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.101343684     -0.004565104      0.001145824
moran.test(res, cont.listw)
## 
##  Moran I test under randomisation
## 
## data:  res  
## weights: cont.listw    
## 
## Moran I statistic standard deviate = 3.0649, p-value = 0.001089
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.101343684      -0.002638522       0.001151024

Test Morana sprawdza, czy wartości reszt są przestrzennie skorelowane, czyli czy sąsiednie powiaty mają podobne wartości reszt. Statystyka testu wynosi 0,1013, a p-value 0.00, co oznacza, że odrzucamy hipotezę zerową, która mówi, że reszty są rozmieszczone losowo. Możemy zatem przyjąć, że reszty w modelu wykazują istotną przestrzenną autokorelację, przez co mogą tworzyć lokalne skupiska wyższych i niższych wartości.

reszty<-factor(cut(res, breaks=c(-100, 0, 100), labels=c("ujemne","dodatnie")))
reszty<-factor(cut(res, breaks=c(-100, 0, 100), labels=c("negative","positive")))
joincount.test(reszty, cont.listw)
## 
##  Join count test under nonfree sampling
## 
## data:  reszty 
## weights: cont.listw 
## 
## Std. deviate for negative = 1.2382, p-value = 0.1078
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic           Expectation              Variance 
##             61.704690             59.013193              4.725001 
## 
## 
##  Join count test under nonfree sampling
## 
## data:  reszty 
## weights: cont.listw 
## 
## Std. deviate for positive = 0.98994, p-value = 0.1611
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic           Expectation              Variance 
##             38.972222             37.013193              3.916184

Nastepnie przeprowadzono Join Count Test, który analizuje ujemne i dodatnie reszty. W przypadku reszt negatywnych statystyka wynosi 1,23, a p-value =0,1078, co wskazuje na brak odrzucenia hipotezy zerowej, dlatego brak jest istotnych klastrów ujemnych. Podobnie jest w przypadku reszt pozytywnych, gdzie statystyka testowa wynosi 0,98, a p-value 0,16, czyli również w tym przypadku brak jest istotnych klastrów dodatnich reszt.

Podsumowując, brak jest silnych skupisk regionów o podobnych resztach, jednakże test Morana pokazuje, że model MNK nie jest w pełni poprawny względem przestrzennym - istnieją istotne zależności przestrzenne, które w tym modelu są niewyjaśniane. Tym samym rozważanie modeli regresji przestrzennej jest jak najbardziej uzasadnione.

Model GWR

bw<-ggwr.sel(eq, data=DaneFinal, coords=crds, family=poisson(), longlat=TRUE)
## Bandwidth: 340.6084 CV score: 29.86086 
## Bandwidth: 550.5657 CV score: 29.9657 
## Bandwidth: 210.8476 CV score: 29.71058 
## Bandwidth: 130.6511 CV score: 29.99905 
## Bandwidth: 260.4118 CV score: 29.77358 
## Bandwidth: 180.2153 CV score: 29.70099 
## Bandwidth: 182.4349 CV score: 29.69951 
## Bandwidth: 190.9971 CV score: 29.69788 
## Bandwidth: 198.5793 CV score: 29.70065 
## Bandwidth: 189.4796 CV score: 29.69776 
## Bandwidth: 189.1548 CV score: 29.69775 
## Bandwidth: 189.0423 CV score: 29.69775 
## Bandwidth: 189.0499 CV score: 29.69775 
## Bandwidth: 189.0502 CV score: 29.69775 
## Bandwidth: 189.0502 CV score: 29.69775 
## Bandwidth: 189.0502 CV score: 29.69775 
## Bandwidth: 189.0502 CV score: 29.69775
model.ggwr<-ggwr(eq, data=DaneFinal, coords=crds, family=poisson(), longlat=TRUE, bandwidth=bw)
model.ggwr
## Call:
## ggwr(formula = eq, data = DaneFinal, coords = crds, bandwidth = bw, 
##     family = poisson(), longlat = TRUE)
## Kernel function: gwr.Gauss 
## Fixed bandwidth: 189.0502 
## Summary of GWR coefficient estimates at data points:
##                              Min.     1st Qu.      Median     3rd Qu.
## X.Intercept.           0.12142705  0.72157172  0.99416684  1.14502860
## Feminizacja           -0.01091336 -0.00576591 -0.00515288 -0.00396955
## log_Imigranci          0.01288642  0.01552752  0.01672913  0.01929564
## wykrywalnosc           0.00233155  0.00391101  0.00485618  0.00577388
## log_rozwody_na10000    0.04813293  0.07688579  0.09275396  0.11219309
## bezrobocie            -0.01056213 -0.00965101 -0.00765726 -0.00530142
## wynagrodzenie_relacja -0.00116723  0.00027860  0.00045897  0.00058340
## log_wydatki_bezp      -0.00262884  0.00945814  0.01102111  0.01246063
## urban                 -0.00403326  0.02845219  0.05591592  0.10180781
## log_zalud             -0.01047042  0.00409560  0.01084745  0.01651792
## log_gosp_jedno         0.16250361  0.18590036  0.20842730  0.24115556
##                              Max.  Global
## X.Intercept.           1.53598371  0.9388
## Feminizacja            0.00046866 -0.0047
## log_Imigranci          0.03474528  0.0192
## wykrywalnosc           0.00684662  0.0044
## log_rozwody_na10000    0.14271023  0.0866
## bezrobocie             0.00125694 -0.0056
## wynagrodzenie_relacja  0.00152739  0.0004
## log_wydatki_bezp       0.01467196  0.0091
## urban                  0.28434515  0.0973
## log_zalud              0.02441462  0.0088
## log_gosp_jedno         0.29470428  0.2013

W modelu pasmo konwergowało w 17 iteracjach od wartości początkowej 340.6 do 189.05. Większość (6) współczynników osiąga znaczący poziom zmienności - przyjmują zarówno wartości dodatnie, jak i ujemne. Wpływ imigracji na przestępczość jest dodatni, co oznacza, że w wiekszości miejsc, wyższa liczba imigrantów powiązana jest z wyższą liczbą czynów nielegalnych. Również wykrywalność, rozwody na 10 000 mieszkańców i odsetek gospodarstw domowych wykazują dodatnie zależności. Marginalny wpływ, bliski zeru obserwowany jest wśród zmiennych, jak: wynagrodzenie_relacja, feminizacja czy bezrobocie.

POW$GWR.log_wydatki_bezp <- model.ggwr$SDF$log_wydatki_bezp
ggplot()+geom_sf(data=POW, aes(fill=GWR.log_wydatki_bezp))

Logarytm wydatków na bezpieczeństwo ma zróżnicowany przestrzennie wpływ na przestępczość w powiatach. W niektórych, jak w przypadku południowo-wschodniej Polski oddziaływanie te jest negatywne, gdzie wzrost wydatków na bezpieczeństwo prowadzi do mniejszej ilości liczby przestępstw. Natomiast w większości powiatów, czyli w centralnej i północnej Polsce wykazywany jest pozytywny wpływ wydatków. Jednakże zmienna ta w skali kraju przyjmuje oszacowanie 0,0091, czyli jej wpływ jest blisku zeru. Dopiero wzrosty wydatków o 1000% przyczyniałyby się do znaczącej zmiany występowania przestępstw.

Współczynnik przy feminizacji w skali kraju ma ujemny wpływ, co może sugerować, że powiaty cechujące się wyższym udziałem kobiet w społeczeństwie mają nizśze poziomy przestępczości. Najsilniejsze negatywne oddziaływanie widoczne jest w północno-zachodniej i południowo-wschodniej części Polski. Efekt jest zróżnicowany przestrzennie, dlatego przyjmuje również niewielkie wartości dodatnie, które obserwowane są głównie w centralnej i północno-wschodniej Polsce. Jednakże maksymalne oszacowanie, to 0.0005, co wskazuje na praktycznie neutarlny wpływ.

Imigranci zwiększają natomiast wartość przestępczości, czyli w powiatach z wyższą liczbą imigrantów popełnianie przestępstw jest częstszym zjawiskiem. Najsilniejszy wpływ widoczny jest w północno-zachodniej Polsce, a najsłabszy w centralnej części kraju, gdzie liczba imigrantów jest największa, ze względu na występowanie dużych aglomeracji, jak Warszawa czy też Łódź, które są atrakcyjnym miejscem do zamieszkania. W północno-zachodniej Polsce liczba przebywających imigrantów jest dosyć niska, nie przekraczając w większości powiatów liczby 1000 imigrantów, jednakże wpływ na przestępczość jest tam największy.

Im więcej rozwodów, tym wyższa liczba przestępstw. Powiaty o większej liczbie rozwodów mogą mieć wyższy poziom przestępczości. Zmienna dot. rozwodów przyjmuje tylko dodatnie wartości, a najsilniejszy efekt widoczny jest w północnej części Polski. Najniższe oddziaływanie obserwowane jest w południowych częściach Polski oraz w powiatach mieszczących się przy granicach kraju. Być może rozpad rodzin powoduje trudniejszą sytuację społczeno-ekonomiczną, co wpływa na wskaźnik przestępczości. Jak i również wytłumaczeniem może być powiązanie rozwodów z występowaniem przestępstw przeciwko rodzinie i opiece, które to również zawierają się w liczbie ogólnie stwierdzonych przestępstw przez Policję.

Urbanizacja ma pozytywny wpływ na zmienną zależną, najsilniejszy efekt widoczny jest w południowo-wschodniej Polsce, a najsłabszy na południowo-zachodniej, w okolicach województwa dolnośląskiego. Minimalny efekt powiązany z ciemnoniebiskim kolorem ma praktycznie neutralny wpływ: -0,004. Województwo podkarpackie oraz lubelskie, które zostały wyszczególnione jasnoniebieskim kolorem należą do wysoko zurbanizowanych terenów i tym samym siła oddziaływania urbanizacji na przestępstwa jest tutaj najwyższa.

Wyższa gęstość zaludnienia zwiększa zmienną zależną, ale efekt ten nie jest jednolity. W południowych regionach wpływ ten jest słabszy, a nawet odwrotny, gdzie minimum osiąga -0.01. Najsilniejszy wpływ widoczny jest w północnej części Polski, gdzie maksymalny efekt wynosi 0.02. Miejsca o większej gęstości zaludnienia charakteryzują się wyższym poziomem przestępczości, ale efekt ten jest nierównomierny przestrzennie. Wyższa gęstość zaludnienia we wschodniej części kraju oddziałuje słabiej na przestępczość niż wyższa gęstość w zachodniej części Polski.

Logarytm odsetka gospodarstw jednoosobowych ma silny pozytywny wpływ na liczbę popełnianych przestępstw. Największe oddziaływanie występuje w południowo-zachodniej części Polski, a dokładniej w województwie dolnośląskim, które cechuje się wysoko zurbanizowanymi i zaludnionymi obszarami. Najsłabszy efekt, jednakże nadal wysoko pozytywny widoczny jest we wschodnich rejonach Polski, w obszarach przygranicznych. Potwierdzają się zatem przesłanki z modelu OLS, gdzie odsetek jednoosobowych gospodarstw domowych powodował wysoki pozytywny efekt na ilość stwierdzanych przestępstw.

Stopa bezrobocia ma zróżnicowany przestrzennie wpływ na liczbę przestępstw, jednakże przeważają wartości ujemne. Oznacza to, że w większości regionów wyższe bezrobocie wiąże się ze spadkiem przestępczości. Najsilniejszy negatywny efekt widoczny jest w południowo-zachodniej i centralnej Polsce, gdzie wraz ze wzrostem stopy bezrobocia ilość popełnianych przestępstw wzrasta. Wpływ dodatni, a bardziej neutralny ze względu na oszacowanie bliskie zeru mają obserwują powiaty znajdujące się w północno-zachodnich terenach kraju.

Wpływ wynagrodzenia na przestępczość, podobnie jak w przypadku bezrobocia jest niewielki. W głównej mierze przeważa jego pozytywne oddziaływanie, które widoczne jest prawie w całej Polsce (wyjątkiem jest północno-zachodnia część kraju, gdzie obserwowane jest negatywne oddziaływanie wynagrodzenia). Dodatnia korelacja jest sprzeczna z intuicją i literaturą przedmiotu, ponieważ zgodnie z teorią przestępczości wyższe wynagrodzenie będzie demotywowało do uciekania się do nielegalnych czynów ze względu na rosnące koszty popełnienia przestępstwa.

Klasteryzacja

# Ile clustrow?
n_clust <- n_clusters(as.data.frame(model.ggwr$SDF[,3:12]),
                      package = c("easystats", "NbClust", "mclust"),
                      standardize = TRUE
)
n_clust
## # Method Agreement Procedure:
## 
## The choice of 4 clusters is supported by 9 (31.03%) methods out of 29 (Silhouette, kl, Hartigan, CCC, Marriot, Tracew, Rubin, Ratkowsky, PtBiserial).
plot(n_clust)

W celu rozszerzenia badania o analizę powiązanych rozkładów postanowiono przeprowadzić klasteryzację. Opieranie się na jednej metodzie poszukiwania klastrów jest mało efektywne, dlatego odwołano się do pakietu NbClust, który przeprowadza wnioskowanie na podstawie 29 metod. Tym samym, aż 9 metod (Silhouette, kl, Hartigan, CCC, Marriot, Tracew, Rubin, Ratkowsky, PtBiserial) wskazuję na wybranie optymalnej liczby klastrów jako 4.

klastry1<-eclust(as.data.frame(model.ggwr$SDF[,3:12]), "kmeans", k=4)

Klastry wyglądają na dobrze podzielone, a metoda K-Means w zadowalający sposób zgrupowała punkty, ponieważ żaden z klastrów nie nachodzi na siebie znacząco.

POW$clust1<-klastry1$cluster
ggplot() + geom_sf(data=POW, aes(fill=clust1))

Wykres przedstawia wyodrębnienie klastrów na podstawie Mapy Polski.

Model błędu przestrzennego

for (i in 1:4) {
  DaneFinal[[paste0("clust", i)]] <- as.integer(klastry1$cluster == i)
}

eq1 = log_przestepstwa~Feminizacja+log_Imigranci+wykrywalnosc+bezrobocie+log_rozwody_na10000+
  bezrobocie+wynagrodzenie_relacja+log_wydatki_bezp+urban+log_zalud+log_gosp_jedno+
  clust2+clust3+clust4

model.sem<-errorsarlm(eq1, data=DaneFinal, cont.listw)
summary(model.sem)
## 
## Call:errorsarlm(formula = eq1, data = DaneFinal, listw = cont.listw)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.54075 -0.17109 -0.02918  0.13660  1.31694 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                         Estimate Std. Error z value  Pr(>|z|)
## (Intercept)            2.4913397  1.0372641  2.4018 0.0163130
## Feminizacja           -0.0147839  0.0088513 -1.6702 0.0948713
## log_Imigranci          0.0519996  0.0181736  2.8613 0.0042195
## wykrywalnosc           0.0134762  0.0019437  6.9334 4.109e-12
## bezrobocie            -0.0126708  0.0121938 -1.0391 0.2987519
## log_rozwody_na10000    0.2787451  0.0741728  3.7581 0.0001712
## wynagrodzenie_relacja  0.0015221  0.0016151  0.9424 0.3459642
## log_wydatki_bezp       0.0234395  0.0134132  1.7475 0.0805521
## urban                  0.2279564  0.1460278  1.5610 0.1185126
## log_zalud              0.0482727  0.0256796  1.8798 0.0601344
## log_gosp_jedno         0.5862674  0.1322082  4.4344 9.232e-06
## clust2                 0.0124953  0.0479337  0.2607 0.7943395
## clust3                 0.0698963  0.0486885  1.4356 0.1511211
## clust4                -0.0014010  0.0527710 -0.0265 0.9788202
## 
## Lambda: 0.21291, LR test value: 8.1095, p-value: 0.0044034
## Asymptotic standard error: 0.071281
##     z-value: 2.9869, p-value: 0.0028178
## Wald statistic: 8.9218, p-value: 0.0028178
## 
## Log likelihood: -38.96349 for error model
## ML residual variance (sigma squared): 0.071255, (sigma: 0.26694)
## Number of observations: 380 
## Number of parameters estimated: 16 
## AIC: 109.93, (AIC for lm: 116.04)

Na bazie klastrowania przypisano obserwacje do klastrów i utworzono zmienne zero-jedynkowe, które uwzględniono w modelu. Tym samym uzyskano kontrole lokalnego podobieństwa, a także autokorelację przestrzenną w składniku losowym. Model wskazuje, że zmienne istotne to zmienne dot. liczby imigrantów, wykrywalności, rozwodów oraz odsetka gospodarstw jednoosobowych. Oddziaływają one silnie na liczbę popełnianych przestępstw, blisko istotności statystycznej znajdują się zmienne dot. feminizacji, wydatków na bezpieczeństwo, czy zaludnienia. W nawiązaniu do modelu OLS zauważalny jest wzrost błędu standardowego wskaźnika urbanizacji. Również niewidoczny jest istotny wpływ klastrów, które oprócz klastra 3 przyjmują wysokie wartości błędów standardowych.

Parametr lambda wskazuje na istotną autokorelacje przestrzenną w składniku losowym, co sugeruje, że uwzględnienie struktury przestrzennej jest konieczne. Wynik testu LR stwierdza, że model SEM jest istotnie lepszy niż klasyczny model MNK.

Podsumowanie

W niniejszej pracy przeprowadzono analizę przestrzenną determinantów przestępczości w Polsce za pomocą regresji ważonej geograficznie (GWR). Model przestrzenny GWR wykazał na istotne zróżnicowanie czynników wpływających na przestępczość w Polsce. Wysoka liczba imigrantów, rozwodów czy odsetka gospodarstw jednoosobowych przyczynia się do zwiększania liczby stwierdzanych przestępstw. Największe oddziaływanie imigracji obserwowano w północno-zachodniej Polsce, mimo stosunkowej niskiej liczby imigrantów w tym regionie. Rozwody najsilniej wpływają na przestępczość w północnej części kraju, a odsetek gospodarstw jednoosobowych ma wyraźnie większy wpływ w regionach Dolnego Śląska. Urbanizacja i zaludnienie również wykazują pozytywne, choć zróżnicowane wpływy. Najsilniejsze oddziaływanie urbanizacji widoczne jest na południowo-wschodnich obszarach Polski, natomiast gęstość zaludnienia ma największe znaczenie w części północnej.

Przeciwny wpływ widoczny w literaturze wykazuje bezrobocie i relacja wynagrodzenia do średniej krajowej. Zgodnie z intuicją wyższa stopa bezrobocia powinna zwiększać skłonność do popełniania nielegalnych czynów, natomiast wyższe wynagrodzenie powinno obniżać intencję do przestępstw. Być może wynika to z charakterystyki mierzenia ogólnej ilości stwierdzonych przestępstw, w których zarówno zawarte są: przestępstwa gospodarcze, przeciwko mieniu, czy w sprawach rodziny i opieki.

Podsumowując, analiza GWR potwierdza istotną przestrzenną zmienność wpływu czynników społeczno-ekonomicznych na przestępczość w Polsce. Tym samym stosowanie strategii prewencyjnych powinno być dostosowywane do specyfiki danych regionów. Wydaję się, że bardziej wartościowe będzie badanie poszczególnych czynów zabronionych, ponieważ inna determinacja może oddziaływać na popełnianie przestępstw gospodarczych, a inna do popełnienia przestępstw przeciwko mieniu.

Literatura

G. Becker, Crime and Pinushment: an economic approach, 1968

G. Dubiec, Wykorzystanie wybranych metod wielowymiarowej analizy porównawczej do oceny poziomu zjawiska przestępczości w powiatach Polski w 2012 roku, 2020

Iwona Bąk, Struktura i typologia przestrzenna przestępczości w Polsce, 2015

J.J. Sztaudynger, M. Sztaudynger, Ekonometryczne modele przestępczości, 2003

K. Kądziołka, Determinanty przestępczości w Polsce. Analiza zależności z wykorzystaniem drzew regresyjnych, 2016

K. Kądziołka, Wpływ wybranych czynników o charakterze społeczno-ekonomicznym na przestępczość przeciwko mieniu w Polsce, 2014

Ł. Gadomski, M. Malaczewski, Zastosowanie dynamicznych modeli panelowych w analizie determinant przestępczości w polsce w latach 2005–2016, 2020

P. Bieniek, S. Cichocki, M. Szczepaniec, Wyniki ekonomiczne a poziom przestępczości - badanie ekonometryczne, 2012

S. Mordwa, Kradzieże w przestrzeni Łodzi, 2011

Tavares JP, Costa AC. Spatial Modeling and Analysis of the Determinants of Property Crime in Portugal. ISPRS International Journal of Geo-Information. 2021; 10(11):731. https://doi.org/10.3390/ijgi10110731

W. Florczak, Makroekonomiczny model przestępczości i systemu egzekucji prawa dla Polski. Specyfikacje równań stochastycznych i rezultaty szacowania parametrów strukturalnych, 2012