ONE-WAY ANOVA analizy wpływuodległości od Paryża na przestępczość we Francji w 1833 roku

PL190120201944

Guerry, Essay on the Moral Statistics of France (1833)
Source: http://datavis.ca/gallery/guerry/guerrydat.html

1. dept wydział Num ID działu Standardowe numery dla działów, z wyjątkiem Korsyki (200).
2. Region Region Char (1) Region Region Francji („N” = „Północ”, „S” = „Południe”, „E” = „Wschód”, „W” = „Zachód”, „C” = „Środek”). Korsyka ma kod „” (brak, NA)
3. Department Departament Char (25) Nazwa departamentu Nazwy departamentów są nazywane zgodnie z użyciem w 1830 roku. Zobacz depts.txt, aby uzyskać alternatywne kodowanie nazw departamentów.
4. Crime_pers Crime_pers Num Muzyka pop. za przestępstwa przeciwko osobom Źródło: A2 (Compte général, 1825-1830)
5. Crime_prop Crime_prop Num Muzyka pop. za przestępstwo przeciwko mieniu Źródło: A2 (Compte général, 1825-1830)
6. Literacy Alfabetyzacja Num Procent odczytu i zapisu Procent poborowych wojskowych, którzy potrafią czytać i pisać Źródło: A2
7. Donations Darowizny Num Darowizny dla biednych Źródło: A2 (Bulletin des lois)
8. Infants Niemowlęta Num Muzyka pop. na nielegalne urodzenie Źródło: A2 (Bureaau des Longitudes, 1817-1821)
9. Suicides Samobójstwa Num Muzyka pop. na samobójstwo Źródło: A2 (Compte général, 1827-1830)
10. MainCity Główne Miasto Char (5) Rozmiar głównego miasta Rozmiar głównego miasta („1: Sm”, „2: Med”, „3: Lg”), stosowany jako surogat dla gęstości zalewania. Duże odnosi się do górnej 10, małe do dolnej 10; wszystkie pozostałe są sklasyfikowane jako średnie. Źródło: A1
11. Wealth Bogactwo Ranga Podatek per capita od nieruchomości osobistych Indeks rankingowy oparty na podatkach od nieruchomości osobistych i ruchomych na mieszkańca Źródło: A1
12. Commerce Handel Ranga Handel i przemysł Handel i przemysł mierzony rangą liczby patentów / populacji Źródło: A1
13. Clergy Kler Ranga Dystrybucja duchowieństwa Dystrybucja duchowieństwa mierzona stopniem liczby księży katolickich w czynnej służbie / ludności Źródło: A1 (Almanach officiel du clergy, 1829)
14. Crime_parents Kryminalni rodzice Ranga Przestępstwa przeciwko rodzicom Przestępstwa przeciwko rodzicom, mierzone stopniem stosunku przestępstw przeciwko rodzicom do wszystkich przestępstw – Średnia z lat 1825–1830 Źródło: A1 (Compte général)
15. Infanticide Dzieciobójca Ranga Liczba dzieciobójstw na jednego mieszkańca Wskaźnik wskaźnikowy liczby dzieciobójstw do liczby ludności – Średnia z lat 1825–1830 Źródło: A1 (Compte général)
16. Donation_clergy Duchowieństwo Ranga Darowizny dla duchowieństwa Wskaźniki rankingowe liczby zapisów i darowizn między żywymi dla ludności – Średnia dla lat 1815–1824 Źródło: A1 (Bull. Des lois, ordunn. D’autorization)
17. Lottery Loteria Ranga Zakład per capita w Royal Lottery Ranking rankingowy wpływów z zakładów loterii królewskiej do liczby ludności — Średnia z lat 1822–1826 Źródło: A1 (Compte rendus par le ministre des finances)
18. Desertion Dezercja Ranga Dezercja wojskowa Desercja wojskowa, stosunek liczby młodych żołnierzy oskarżonych o dezercję do siły kontyngentu wojskowego, minus deficyt spowodowany niewystarczalnością dostępnych polan – średnia z lat 1825–1827. Źródło: A1 (Compte du ministere du guerre, 1829 etat V)
19. Instruction Instrukcja Ranga Rankingi instrukcji zapisane z mapy instrukcji Guerry’ego. Uwaga: jest to odwrotnie związane z umiejętnością czytania i pisania (zgodnie z definicją tutaj).
20. Prostitutes Prostytutki Num Prostytutki w Paryżu Liczba prostytutek zarejestrowanych w Paryżu w latach 1816–1834, sklasyfikowanych według departamentu ich urodzenia Źródło: Parent-Duchatelet (1836), De la prostitution en Paris
21. Distance Dystans Num Odległość do Paryża (km) Odległość każdego centroidu departamentu do centroidu Sekwany (Paryż) Źródło: obliczone na podstawie centroidów departamentów
22. Area Powierzchnia Num Obszar (1000 km^2) Źródło: Angeville (1836)
23. Pop1831 Pop1831 Num 1831 populacja ludności w 1831 roku, pochodzi z Angeville (1836), Essai sur la Statistique de la Populacja français w 1000s.

In [2]:
import pandas as pd

df = pd.read_csv('c:/1/Guerry.csv', index_col='Department')
del df['Unnamed: 0']
df.head(4)
Out[2]:
dept Region Crime_pers Crime_prop Literacy Donations Infants Suicides MainCity Wealth Crime_parents Infanticide Donation_clergy Lottery Desertion Instruction Prostitutes Distance Area Pop1831
Department
Ain 1 E 28870 15890 37 5098 33120 35039 2:Med 73 71 60 69 41 55 46 13 218.372 5762 346.03
Aisne 2 N 26226 5521 51 8901 14572 12831 2:Med 22 4 82 36 38 82 24 327 65.945 7369 513.00
Allier 3 C 26747 7925 13 10973 17044 114121 2:Med 61 46 42 76 66 16 85 34 161.927 7340 298.26
Basses-Alpes 4 E 12935 7289 46 2733 23018 14238 1:Sm 76 70 12 37 80 32 29 2 351.399 6925 155.90

4 rows × 22 columns

Celem naszego badanie jest sprawdzenie czy istnieje istotny statystycznie wpływ miasta Paryża na poziom przestępczości przeciwko osobom (Crime_pers). Podzielę wszystkie departamenty na 4 strefy odległości od miasta Paryża (Distance) i porówna poziom przestępczości przeciwko osobom.

– h0: poziom przestępczości (średnia) przeciwko osobom jest taki sam pomiędzy strefami
– h1: poziom przestępczości (średnia) przeciwko osobom jest statycznie inny pomiędzy strefami

In [9]:
df['Class_of_Distance']=pd.qcut(df['Distance'], 4,labels=["do 120 km","120-200 km","200-280 km","280-539 km"])
In [12]:
pd.pivot_table(df, index= ['Class_of_Distance'], values= "Distance", aggfunc= ['min','max', 'count'])
Out[12]:
min max count
Distance Distance Distance
Class_of_Distance
do 120 km 0.000 119.718 22
120-200 km 126.378 199.167 21
200-280 km 202.065 283.810 21
280-539 km 291.624 539.213 22

Nazwałęm te klasy dopiero po tym jak komputer je wskazał. Innymi słowy qcut – podzielil miasteczka na 4 klasy po równo, ja dodałem odpowiednie nazwy klas. Mając klasy Tworzę odpowiednią data frame aby utworzyć One-way ANOVA.

In [14]:
df[['Crime_pers','Class_of_Distance']].head(3)
Out[14]:
Crime_pers Class_of_Distance
Department
Ain 28870 200-280 km
Aisne 26226 do 120 km
Allier 26747 120-200 km

Test ANOVA

In [16]:
import statsmodels.api as sm
from statsmodels.formula.api import ols

model_H = ols('Crime_pers ~ C(Class_of_Distance)', data=df).fit()

anova_table = sm.stats.anova_lm(model_H, typ=2)
print(anova_table)
                            sum_sq    df          F    PR(>F)
C(Class_of_Distance)  1.408371e+09   3.0  11.392976  0.000003
Residual              3.378877e+09  82.0        NaN       NaN

Interpretacja: Wartość P uzyskana z analizy ANOVA jest znacząca (P

ANOVA pokazałe, że istnieje statystycznie istotna różnice w poziomie przestępczości. Niestety ANOVA nie wskazuje, które grupy obszary istotnie różnią się od siebie.
Aby poznać pary obszarów różniące się istotnie poziomem prtzestępczości, należy przeprowadzić analizę wielokrotnego porównania par (porównanie post-hoc) za pomocą testu HSD Tukeya.

Testu HSD Tukeya

In [18]:
from statsmodels.stats.multicomp import pairwise_tukeyhsd

m_comp = pairwise_tukeyhsd(endog=df['Crime_pers'], groups=df['Class_of_Distance'], alpha=0.05)
print(m_comp)
       Multiple Comparison of Means - Tukey HSD,FWER=0.05      
===============================================================
  group1     group2     meandiff     lower      upper    reject
---------------------------------------------------------------
120-200 km 200-280 km   -6203.0   -11398.3884 -1007.6116  True 
120-200 km 280-539 km -10450.2446 -15586.2551 -5314.234   True 
120-200 km do 120 km   -1923.1537  -7059.1642 3212.8569  False 
200-280 km 280-539 km  -4247.2446  -9383.2551  888.766   False 
200-280 km do 120 km   4279.8463   -856.1642  9415.8569  False 
280-539 km do 120 km   8527.0909   3451.1527  13603.0291  True 
---------------------------------------------------------------
In [37]:
m_comp.plot_simultaneous()    # Plot group confidence intervals
plt.vlines(x=20000,ymin=-0.5,ymax=3.5, color="red", alpha=0.8, linestyle='--')
Out[37]:
<matplotlib.collections.LineCollection at 0x14963a6ecc0>

Test Tukey HSD wykazał statystycznie istotną różnicę w przestępczości pomiędzy trzeba parami obszarów0:

– “120-200 km” i “200-280 km”
– “120-200 km” i “280-590 km”
– “280-590 km” i “do 120 km”

Powyższe wyniki z Tukey HSD sugerują, że oprócz wymienionych zestawów, wszystkie inne porównania par zachowują hipotezę zerową czyli w pozostałych obszarach nie istnieją istotne różnice statystyczne w poziomie przestępczości przeciwko osobom.

Sprawdzenie spełnienia warunków ANOVA

Warunki:

  1. wartości rezydualne mają rozkład normalny (test Shapiro Wilksa)
  2. wariancje w grupach są jednorodne (test Levene lub Bartlett)
  3. obserwacje są prowadzone niezależnie od siebie

Test Levene’a Sprawdzenie jednorodności wariancji

Hipoteza zerowa : grupy z populacji mają równe wariancje.

In [20]:
PKS = pd.pivot_table(df, index = 'Department', columns = 'Class_of_Distance', values='Crime_pers')
PKS.head(4)
Out[20]:
Class_of_Distance do 120 km 120-200 km 200-280 km 280-539 km
Department
Ain NaN NaN 28870.0 NaN
Aisne 26226.0 NaN NaN NaN
Allier NaN 26747.0 NaN NaN
Ardeche NaN NaN 9474.0 NaN
Test Levene’a nie toleruje pustych wartości NaN

In [21]:
P01=PKS['do 120 km'].dropna(how='any')
P02=PKS['120-200 km'].dropna(how='any')
P03=PKS['200-280 km'].dropna(how='any')
P04=PKS['280-539 km'].dropna(how='any')

Wykresy przedstawiają statystyki poziomu przestępczości w zadanych czterech obszarach.

In [29]:
PKS.plot.kde()
Out[29]:
<matplotlib.axes._subplots.AxesSubplot at 0x149634c7da0>
In [31]:
import matplotlib.pyplot as plt
PKS.boxplot(column=['do 120 km', '120-200 km', '200-280 km', '280-539 km'], grid=False)
Out[31]:
<matplotlib.axes._subplots.AxesSubplot at 0x149637838d0>
In [30]:
import scipy.stats as stats
w,p = stats.levene(P01,P02,P03,P04)
print("Value:   ",w)
print("p-value: ",p)
Value:    1.2094373013491309
p-value:  0.31156953536619264

Ponieważ wartość P nie jest nieznacząca (p>0.05), nie odrzucamy hipotezy zerowej – czyli grupy wiekowe mają jednorodne wariancje.

Test Shapiro-Wilk Sprawdzenia normalności rozkładu reszt

Hipoteza zerowa : wartości rezydualne mają rozkład normalny.

In [25]:
import scipy.stats as stats
import numpy as np


w, p = stats.shapiro(model_H.resid)
print("Value:   ",w)
print("p-value: ",np.round(p, decimals=2))
Value:    0.9815720915794373
p-value:  0.26

Ponieważ wartość P nie jest znaczące bo większa od współczynnika ufności 0.05 (p>0.05), nie ma podstaw do odrzucenia hipotezy zerowej – wartości rezydualne mają rozkład normalny.

Przyjrzymy się bliżej resztą modelu.

In [26]:
from statsmodels.graphics.gofplots import qqplot
from matplotlib import pyplot

qqplot(model_H.resid, line='s')
pyplot.show()

Faktynie reszty zachowują prawie idealny rozkłąd normalny.

 PL080120201039