Co ma wpływ na wysokość zarobków?¶
1. Wstęp¶
Chciałbym dzisiaj poruszyć nieco nietypowy temat i spróbować odpowiedzieć na pytanie: co ma wpły na wysokość zarobków? Wybór tematu nie jest przypadkowy. Pozwoli mi on ukazać wady i zalety, a także pułapki wnioskowania statystycznego opartego o model regresji liniowej. Zaczynamy!
Na wstępie kilka uwag. Wpis ten znacząco różni się od poprzednich projektów opublikowanych na blogu z trzech powodów:
- Skupia się on w nieco większym stopniu na wnioskowaniu parametrycznym w oparciu o model, niż na "szlifowaniu" wyniku.
- Jest to pierwszy projekt na blogu napisany od początku* do końca w Pythonie.
- Testuję nowy sposób na publikację projektów - notatniki Jupyter. Do tej pory pisałem kod w IDE (R Studio, lub Spyder), dorabiałem do tego opis w MS Word bazując na kometarzam z kodu, a następnie składałem wszystko w jedną całość w edytorze wpisów WordPress. Praca z notebookami Jupyter wygląda zupełnie inaczej. Wszystko znajduje się w jednym miejscu. Komentarze i kod wraz z wynikami są zapisywane w jednym notatniku, który jest gotowy do publikacji do formatu HTML. Opublikowany kod HTML wystarczy skopiować do WordPress i gotowe. Pozwala to zaoszczedzić czas, którego mi ostatnio bardzo brakuje. Ciekaw jestem czy taka forma prezentacji przypadnie Wam do gustu :)
* Przy wczytywaniu danych wspomogłem się językiem R. Udostępnione dane z których korzystałem były zapisane w formacie ".dta", które już kiedyś wczytywałem w R-ce.
1.1. Opis problemu¶
Na świecie istnieją różnice w zarobkach i chyba nikogo nie trzeba do tego przekonywać. W różnych regionach świata, różni ludzie, o różnym wykształceniu otrzymują różne stawki za wykonywaną pracę. Czy zastanawiałeś się kiedyś co w największym stopniu wpływa na wysokość zarobków? Czy jest to zdybyte doświadczenie (wiek danej osoby)? A być może coś mniej oczywistego, np. płeć i rasa? No i w zasadzie jaki wpływ na pensje, niezależnie od płci i pochodzenia ma edukacja? W tym wpisie postaram się odpowiedzieć na te i kilka innych pytań dotyczących zarobków :)
1.2. Opis zbioru danych¶
By odpowiedzieć wszystkie pytania wymienione w poprzednim punkcie potrzebna jest dogłębna analiza i wnioskowanie parametryczne. Analiza ma sens jedynie wtedy, gdy dysponujemy odpowiednio wielką próbką danych o wysokiej jakości. Pozyskanie jej było więc dla mnie priorytetem.
Źródło zbioru
Po rozważeniu kilku opcji zdjecydowałem się użyć zbioru, który był opisywany w znakomitej książce autorstwa Andrew Gelmana i Jennifer Hill "Data Analysis Using Regression and Multilevel/Hierarchical Models".
Pełen zbiór jest dostępny na stronie wydawnictwa, pod linkiem: klik.
Opis zbioru
Zbiór powstał w oparciu o wyniki ankiet przeprowadzonych na mieszkańcach USA w latach 90-tych ubiegłego wieku. Ich celem było badanie wpływu poszczególnych czynników społecznych i środowiskowych na poziom życia mieszkańców Stanów Zjednoczonych.
Przedstawiony w zbiorze poziom zarobków ankietowanego oznacza jego roczny przychód, wyrażony w dolarach. Jeśli chcielibyśmy porównać wyniki widoczne w danych do aktualnej sytuacji gospodarczej, to musielibyśmy wziąć pod uwagę m.in. inflację. Na potrzeby tego projektu nie jest to konieczne i zostawię te wartości niezmienione. Bardziej zależy mi na porównaniu wpływu poszczególnych czynników na wysokość zarobków (w momencie przeprowadzenia ankiety), niż na odnoszeniu się do aktualnej sytuacji w USA.
Warto zaznaczyć, że uzyskane wyniki mogą, ale nie muszą być miarodajne dla innych regionów świata, w tym dla Polski. Gdybyśmy przeprowadzili podobne ankiety w naszym kraju i zinterpretowali ich wyniki, mogłoby sie okazać, że wyniki analiz są zupełnie różne. Zastrzegam zatem, że wszystkie wnioski zawarte w tym wpisie dotyczyć będą jedynie badanej populacji USA :)
2. Przygotowanie zbioru¶
Jak wczytać dane?
Ten projekt bazuje na danych zawartych w pliku "heights.dta". Jeśli chcielibyście wczytaj dane w R, to możecie to zrobić dopiero po instalajci biblioteki "readstata13" (której również ja użyłem):
install.packages("readstata13")
require("readstata13")
dat <- read.dta13("heights.dta")
Tak wczytany zbiór wtępnie obrobiłem, usuwając z niego obserwacje zawierające brakujące wartości. Dla purystów dodam, że Pythonowy Pandas również umożliwia czytanie plików dta. Do dalszej analizy będę używać już tylko i wyłącznie Pythona :)
2.1. Wczytanie podstawowych bibliotek i pliku z danymi¶
Wczytuję niezbędne biblioteki:¶
import pandas as pd # Biblioteka do analizy danych - dataframe-y w Python-ie
import numpy as np # Zestaw narzędzi do obliczeń numerycznych
import matplotlib.pyplot as plt # Wizualizacja danych - podstawowa biblioteka
import seaborn as sns # Wizualizacja danych - nakładka upiększająca :)
import statsmodels.formula.api as sm # Modele statystyczne, w tym regresja liniowa
import scipy.stats as stats # Biblioteka niezbędna m.in. do przeprowadzania testów statystycznych
Zmieniam sposób wyświetlania danych numerycznych:¶
Precyzja do dwóch miejsc po przecinku będzie wystarczająca. Oczywiście dotyczy ona jedynie "drukowanych" wartości.
pd.set_option('float_format', '{:.2f}'.format)
Wczytuję plik i zmieniam nazwy kolumn:¶
df = pd.read_csv('data/zarobki.txt')
df.columns = ['zarobki','wzrost','plec','rasa','lata_edukacji','wiek']
2.2. Zmiana typów zmiennych i jednostek¶
df.dtypes
df.rasa.unique()
Jeśli chodzi o typy zmiennych to wyglądają one następująco:¶
- 'zarobki', 'wzrost', 'lata_edukacji' i 'wiek' - zmienne numeryczne,
- 'plec' i 'rasa' - zmienne kategoryczne (Python interpretuje je jako "object").
Już na pierwszy rzut oka widać, że warto wprowadzić kilka drobnych korekt:¶
- Zaokrąglić wartości zmiennej 'zarobki' i zmienić jej jednostkę na tys. dolarów (dzięki mniejszym liczbom, wyniki będą bardziej przejrzyste).
- Zmienić jednostkę zmiennej 'wzrost' z cali na cm - będzie to pomocne przy interpretacji modelu.
- Zaokrąglić wartości zmiennej 'wzrost', do pełnego cm.
- Zmienić wartości jakie przyjmują zmienne: 'plec' i 'rasa'.
- Zmienić typ zmiennych: 'plec', 'rasa', na kategoryczne.
# 1
df.zarobki = (df.zarobki/1000).round(3)
# 2,3
df.wzrost = (df.wzrost * 2.54).round(0).astype(int)
# 4
df.plec.replace(['male','female'],['mężczyzna','kobieta'], inplace=True)
df.rasa.replace(['white','other','hispanic','black'],['biała','inna','hiszpańska','czarna'], inplace=True)
df.rasa = df.rasa.astype('category')
df.plec = df.plec.astype('category')
Podgląd pierwszych 10 obserwacji ze zbioru:¶
df.head(10)
Weryfikacja zmiennych kategorycznych¶
Sprawdźmy, czy nie ma żadnych niespodzianek w zmiennych kategorycznych (np. błędnie wpisana nazwa kategorii).
print('Płeć - kategorie:')
df.plec.unique()
print('Rasa - kategorie:')
df.rasa.unique()
print('Typy zmiennych:')
df.dtypes
Ok. Teraz całośc wygląda znacznie lepiej :)
3. Eksploracyjna analiza danych¶
3.1. Sprawdzenie podstawowych statystyk w celu wykrycia braków w danych.¶
3.1.1. Statystyki zmiennych numerycznych:¶
df.describe()
3.1.2. Analiza braków danych:¶
summary = pd.DataFrame(df.dtypes, columns=['Dtype'])
summary['Nulls'] = pd.DataFrame(df.isnull().any())
summary['Sum_of_nulls'] = pd.DataFrame(df.isnull().sum())
summary['Per_of_nulls'] = round((df.apply(pd.isnull).mean()*100),2)
summary.sort_index(inplace=True)
summary
Żadnych braków w danych, 1379 obserwacji. Jest nieźle :)
3.2. Dodanie nowej zmiennej - 'wykształcenie'¶
Pierwszą rzeczą która przyszła mi do głowy w poprzednim kroku było dodanie nowej zmiennej, która określałaby poziom edukacji. Coś na wzór polskiego poziomu wykształcenia (podstawowe, średnie, wyższe). Taka zmienna mogłaby wnieść coś nowego do analizy. Punktem wyjścia do podziału zmiennej 'lata_edukacji' na koszyki, powinny być wartości graniczne koszyków. By dokonać ich estymacji muszę się bliżej przyjrzeć rozkładowi zmiennej 'lata_edukacji' i dowiedzieć się nieco więcej o amerykańskim systemie kształcenia z Wikipedii :)
plt.figure(figsize=(15,7))
sns.set(font_scale=1.5)
sns.kdeplot(df.lata_edukacji, shade=True).set(title="Wykres gestosci - 'lata_edukacji'")
plt.show()
plt.figure(figsize=(15,7))
sns.set(font_scale=1.5)
sns.distplot(df.lata_edukacji, bins=20, kde=False).set(title="Histogram - 'lata_edukacji'", ylabel ="Liczba osob")
plt.show()
Widać że zmienna dzieli się na 3 lub 4 'grubsze koszyki'. Zakładam, że to przyczyna systemu edukacyjnego funkcjonującego w USA. Zgodnie z tamtejszym podziałem poziomy edukacji można podzielić następująco:
- szkoła podstawowa: 6-8 lat (średnio 7),
- szkoła średnia: 4-6 lat (średnio 5, sumarycznie 12),
- szkoła policealna (opcjonalnie): 2 lata (sumarycznie 14),
- szkoła wyższa: 4-6 lat (średnio 5, sumarycznie 19).
Powyższe założenia mniej więcej pokrywają się z tym co zostało zaprezentowane na histogramie. Dodaję zatem do zbioru dodatkową zmienną.
df['wyksztalcenie'] = pd.cut(df.lata_edukacji, bins=[0,8,12,14,19], labels=['szkoła podstawowa', 'szkoła średnia', 'szkoła policealna', 'szkoła wyższa'])
plt.figure(figsize=(15,7))
sns.set(font_scale=1.5)
sns.countplot(df.wyksztalcenie, order=df.wyksztalcenie.value_counts().index, palette='Blues').set(xlabel="", ylabel="Liczba osob", title="Liczba osób z danym wykształceniem")
plt.show()
Liczba osób w poszczególnych kategoriach:
print(df.wyksztalcenie.value_counts())
3.3. Analiza liczebności atrybutów kategorycznych¶
Celem tego kroku jest ustalenie liczby obserwacji w wybranych kategoriach. Informacja ta może okazać się pomocna przy dalszej analizie, np. przy badaniu zależności pomiędzy atrybutami.
3.3.1. Płeć¶
plec = pd.DataFrame(df.plec.value_counts()).copy()
plec['Percent'] = ((df.plec.value_counts()/df.shape[0])*100).round(2)
plec
Pierwsze zaskoczenie. Mamy około 65% więcej kobiet niż mężczyzn w naszym zbiorze.
plt.figure(figsize=(15,7))
sns.set(font_scale=1.5)
sns.countplot('plec',data=df, palette='Blues').set(xlabel="", ylabel='Liczba obserwacji', title="Liczba kobiet vs liczba mężczyzn", xticklabels=['Kobiety', 'Mężczyzni'])
plt.show()
3.3.2. Rasa¶
rasa = pd.DataFrame(df.rasa.value_counts()).copy()
rasa['Percent'] = ((df.rasa.value_counts()/df.shape[0])*100).round(2)
rasa
Spoglądając ma powyższe proporce mam wątpliwości co do tego, czy próba odzwierciedla rzeczywiste proporcje ludności zamieszkującej USA :)
plt.figure(figsize=(15,7))
sns.set(font_scale=1.5)
sns.countplot('rasa',data=df, palette='Blues').set(xlabel="", ylabel='Liczba obserwacji', title="Liczba osób z poszczególnych ras w zbiorze")
plt.show()
3.4. Analiza atrybutów numerycznych¶
Prócz wyznaczenia podstawowych statystyk dotyczących zmiennych numerycznych, przeprowadzę też analizę ich rozkładów. W kolejnych krokach będzie ona pomocna z dwóch powodów:
Niektóre metody mają założenia dotyczące rozkładu. Przykładowo, podczas badania korelacji musimy pamiętać, że niektóre metody wyznaczania korelacji zakładają normalność rozkładów (np. współczynnik korelacji Pearsona).
Dzięki weryfikacji rozkładów będę mógł zweryfikować czy istnieją jakiekolwiek wartości odstające dla danej zmiennej.
3.4.1. Zarobki¶
Podstawowe statystyki:¶
summary = pd.DataFrame(df.zarobki.describe().transpose())
summary = summary.transpose()
summary['median'] = df.zarobki.median()
summary
Pierwsze wnioski:
- Średnia jest znacznie większa niż mediana. Może to świadczyć o prawoskośności.
- Wartość minimalna rocznych zarobków, to ok. -$100. Nie wiem z czego to wynika ale dla pewności sprawdzę ile jest obserwacji ujemnych. Jeśli skala jest niewielka, to potraktuję je jako wartości odstające i usunę je.
print("Liczba obserwacji odstających: " + str(df[df.zarobki<0].shape[0]))
df[df.zarobki<0]
Mamy więc raptem 11 takich obserwacji. Są to głównie mężczyźni, w większości w wieku 22-25 lat. Studenci? Być może. Nie mniej, dla uporządkowania danych, pozbędę się ich w punkcie 4.
Badanie rozkładu:¶
plt.figure(figsize=(15,7))
sns.set(font_scale=1.5)
sns.distplot(df.zarobki, bins=30, kde=False).set(title="Histogram - 'zarobki'", ylabel ="Liczba osób")
plt.show()
Obserwacje rozkładają się mocno niesymetrycznie, względem wartości średniej. Rozkład jest prawoskośny. W praktyce oznacza to, że w zbiorze znajduję się więcej osób, które zarabiają poniżej średniej. Średnia jest zawyżana przez obserwacje odsatjące (coś mi mówi, że podobnie wygląda sytuacja zarobków w Polsce). W punkcie 4 będę chciał się ich pozbyć, gdyż regresja liniowa jest wrażliwa na obserwacje odstające. Do ustalenia liczby wartości odstających użyję wykresu pudełkowego i metody "1.5 rozstępu międzykwartylowego".
Warto również zaznaczyć, że powyższy rozkład nie jest rozkładem normalnym. Podczas badania korelacji pomiędzy zmienną 'zarobki' a pozostałymi zmiennymi numerycznymi, nie powinienem używać współczynnika korelacji Pearsona.
plt.figure(figsize=(15,7))
sns.set(font_scale=1.5)
sns.boxplot(df.zarobki, palette='Blues').set(title="Wykre pudełkowy - 'zarobki'")
plt.show()
Wartości odstające:¶
Kropki na powyższym wykresie są outlier-ami. Zgodnie z założeniami wykresu pudełkowego, wartością odstająca jest obserwacja leżąca powyżej 1,5 wartości rozstępu międzykwartylowego (Q3-Q1, obejmuje on 50% procent wszystkich typowych obserwacji) od Q1 (dolne wartości odstające), lub Q3 (górne wartości odstające). Zbiór zawiera więc co najmniej kilkanaście obserwacji, których będę musiał się pozbyć.
Wartość progową, powyżej której obserwacja będzie traktowana jako odstająca, można łatwo obliczyć. Wystarczy użyc powyższej reguły.
q1 = df.zarobki.quantile(0.25)
q3 = df.zarobki.quantile(0.75)
q3+(1.5*(q3-q1))
Zatem wartością graniczną zarobków jest roczny przychód na poziomie $95.458k. Sprawdzę jeszcze ile jest takich obserwacji i czy ich usunięcie nie uszczupli zbytnio zbioru.
population = df.shape[0]
outliers = df[df.zarobki>95.45775].shape[0]
print('Liczba wszystkich obserwacji: ' + str(population))
print('Liczba obserwacji odstających: ' + str(outliers))
print('Odsetek obserwacji odstających w badaniej próbie: ' + str(round(outliers/population*100,2)) + '%')
Nie jest źle. Ubytek na poziomie 3,41% nie powinien negatywnie wpłynąć na dalszą analizę.
3.4.2. Wzrost¶
Podstawowe statystyki:¶
summary = pd.DataFrame(df.wzrost.describe().transpose())
summary = summary.transpose()
summary['median'] = df.wzrost.median()
summary
Coś mi jednak mówi, że w przypadku wzrostu tego typu podsumowanie nie jest miarodajne :) Wzrost jest mocno powiązany z płcią (tego chyba nie muszę dowodzić, więc pozwolę sobie na takie założenie a priori ;-) ), co mocno zaburza wyniki analizy. Wpływa to negatywnie na dalsze wnioski. Dla przykładu: spoglądając na powyższe dane można by stwierdzić, że wykres jest lekko prawoskośny (średnia jest nieco większa od mediany). Zaciera to jednak prawdziwy obraz danych. Wykres jest lekko prawoskośny, ponieważ statystycznie kobiety są nieco niższe od mężczyzn, a w badanej próbie to kobiety stanowią większość.
Na potrzeby dalszej analizy zmiennej "wzrost", oddzielam kobiety od mężczyzn.
wzrost_summary = df.groupby('plec')['wzrost'].describe()
wzrost_summary['median'] = df.groupby('plec')['wzrost'].median().astype(float)
wzrost_summary
Badanie rozkładu:¶
plt.figure(figsize=(15,7))
sns.set(font_scale=1.5)
sns.distplot(df[df.plec=='mężczyzna'].wzrost, bins=25, kde=False, label="Mężczyźni").set(title="Histogram - 'wzrost'", ylabel ="Liczba osob")
sns.distplot(df[df.plec=='kobieta'].wzrost, bins=25, kde=False, label="Kobiety").set(title="Histogram - 'wzrost'", ylabel ="Liczba osob", xlabel='wzrost [cm]')
plt.legend()
plt.show()
W tym przypadku o rozkładach dla obu kategorii niewiele można powiedzieć. Widać oczywiście podstawową różnicę: statystycznie mężczyźni są wyższi od kobiet. Jeśli chodzi o rozkład wzrostu dla kobiet, to nieco przypomina on rozkład normalny. Jest jednak za wcześnie by to potwierdzić.
plt.figure(figsize=(15,7))
sns.set(font_scale=1.5)
sns.boxplot('plec', 'wzrost', data=df, palette='Blues').set(title="Wykre pudełkowy - 'zarobki'")
plt.show()
Wartości odstające:¶
Używając reguły "1.5 wartości rozstępu międzykwartylowego" mogę z łatwością stwierdzić, że zbiór zawiera kilku nietypowo niskich mężczyzn, oraz kilka nietypowo wysokich i niskich kobiet. Rozkład wzrostu dla kobiet nieco przypominał rozkład normalny, dlatego przetestuję na nich inną metodę filtrowania wartości odstających.
Reguła trzech sigm, bo o niej mowa, mówi że 99,7% wartosci danej cechy leży w odległosci trzech sigm (odchyleń standardowych) od wartosci oczekiwanej. Reguła ta ma jednak swoje założenia. Najważniejsze z nich mówi, że dotyczy ona jedynie zmiennych o rozkładzie normalnym. Muszę zatem w pierwszej kolejności sprawdzić rozkład zmiennej 'wzrost' dla obu płci.
Są co najmniej 2 znane mi sposoby na sprawdzenie rozkładu: test normalnosci i analiza histogramu. Histogram znajduje się powyżej i nie jestem w stanie jednoznacznie stwierdzić czy założenie normalności jest spełnione przynajmniej dla kobiet. Poza tym analiza histogramu jest metodą subiektywną. Wolę bazować na faktach, zatem zastosuję test normalności.
Test normalności rozkładu:
Stawiam dwie hipotezy:
- H0 - hipoteza zerowa zakłada, że różnica pomiędzy badanymi parametrami wynosi zero. W tym przypadku będzie to oznaczać, że zmienna ma rozkład normalny.
- H1 - hipoteza alternatywna - przyjmowana jest wtedy, gdy H0 jest odrzucana. Oznaczać będzie, że zmienna nie ma rozkładu normalnego.
Jak działa test?
Jeżeli p_val jest bardzo małe (<alfa), to istnieje małe prawdopodobienstwo popełnienia błędu I rodzaju (odrzucenia hipotezy 0, gdy jest ona prawdziwa).
Jako alfa (poziom istotnosci testu, będący maksymalnym dopuszczalnym prawdopodobieństwem popełnienia błędy I rodzaju, który akceptuję) przyjmuję "klasyczną" wartosc 0.05. Jesli p-value < alfa, to mogę odrzucić H0, gdyż prawdopodobieństwo popełnienia błędu jest akceptowalnie niskie.
print("Test normalności rozkładu zmiennej 'wzrost' dla kobiet:")
if stats.normaltest(df[df.plec=='kobieta'].wzrost)[1]<0.05:
print('\t-Rozkład nie jest rozkładem normalnym')
else:
print('\t-Rozkład jest rozkładem normalnym')
print("Test normalności rozkładu zmiennej 'wzrost' dla mężczyzn:")
if stats.normaltest(df[df.plec=='mężczyzna'].wzrost)[1]<0.05:
print('\t-Rozkład nie jest rozkładem normalnym')
else:
print('\t-Rozkład jest rozkładem normalnym')
Z maksymalnym prawdopodobieństwem popełnienia błędu równym 5%, w obu przypadkach odrzucam hipotezę zerową. Zgodnie z wynikami testu, rozkład wzrostu dla obu płci nie jest rozkładem normalnym. Nie mogę więc w tym przypadku zastosować reguły trzech sigm. Pozostaje mi zatem sprawdzona reguła "1.5 wartości rozstępu międzykwartylowego".
Zgodnie z wykresem pudełkowym dla zmiennej 'wzrost' mamy kilka nietypowych obserwacji:
- kilku niskich mężczyzn
- kilka niskich i wysokich kobiet
Obliczam zatem dolną wartość progową zmiennej 'wzrost' dla mężczyzn:
q1 = df[df.plec=='mężczyzna'].wzrost.quantile(0.25)
q3 = df[df.plec=='mężczyzna'].wzrost.quantile(0.75)
print('Dolna wartość: ' + str(q1-(1.5*(q3-q1))))
Obliczam górą i dolną wartość progową zmiennej 'wzrost' dla kobiet:
q1 = df[df.plec=='kobieta'].wzrost.quantile(0.25)
q3 = df[df.plec=='kobieta'].wzrost.quantile(0.75)
print('Dolna wartość: ' + str(q1-(1.5*(q3-q1))))
print('Górna wartość: ' + str(q3+(1.5*(q3-q1))))
3.4.3. Lata edukacji¶
Podstawowe statystyki:¶
summary = pd.DataFrame(df.lata_edukacji.describe().transpose())
summary = summary.transpose()
summary['median'] = df.lata_edukacji.median()
summary
Badanie rozkładu:¶
# Histogram
plt.figure(figsize=(15,7))
sns.set(font_scale=1.5)
sns.distplot(df.lata_edukacji, bins=10, kde=False).set(title="Histogram - 'lata_edukacji'", ylabel ="Liczba osób")
plt.show()
# Wykres pudełkowy
plt.figure(figsize=(15,7))
sns.set(font_scale=1.5)
sns.boxplot(df.lata_edukacji, palette='Blues').set(title="Wykre pudełkowy - 'lata_edukacji'")
plt.show()
Istnieje kilka osób które swoją edukację skończyły wyjątkowo wcześnie.
Wartości odstające:¶
q1 = df.lata_edukacji.quantile(0.25)
q3 = df.lata_edukacji.quantile(0.75)
print('Dolna wartość: ' + str(q1-(1.5*(q3-q1))))
print('Górna wartość: ' + str(q3+(1.5*(q3-q1))))
Liczba obserwacji odstających:
df[df.lata_edukacji<7.5].shape[0]
3.4.4. Wiek¶
Podstawowe statystyki:¶
summary = pd.DataFrame(df.wiek.describe().transpose())
summary = summary.transpose()
summary['median'] = df.wiek.median()
summary
Badanie rozkładu:¶
plt.figure(figsize=(15,7))
sns.set(font_scale=1.5)
sns.distplot(df.wiek, bins=25, kde=False).set(title="Histogram - 'wiek'", ylabel ="Liczba osób")
plt.show()
plt.figure(figsize=(15,7))
sns.set(font_scale=1.5)
sns.boxplot(df.wiek, palette='Blues').set(title="Wykres pudełkowy - 'wiek'")
plt.show()
Wartości odstające:¶
q1 = df.wiek.quantile(0.25)
q3 = df.wiek.quantile(0.75)
print('Dolna wartość: ' + str(q1-(1.5*(q3-q1))))
print('Górna wartość: ' + str(q3+(1.5*(q3-q1))))
Liczba obserwacji odstających:
df[df.wiek>88].shape[0]
4. Operacje na zbiorze¶
Mając podstawowe informacje o zbiorze, mogę przejść do naniesienia ostatnich korekt.
4.1. Usunięcie wartości odstających¶
Regresja liniowa, którą wykorzystamy do analizy i testowania hipotez jest algorytmem wrażliwym na braki danych i wartości odstające, dlatego tak ważne jest "odsianie" obserwacji nietypowych.
Zgodnie w wnioskami z punktu 3, usuwam ze zbioru osoby o nietypowych zarobkach, wzroście, wieku i liczbie lat edukacji.
4.1.1. Zarobki¶
df = df[(df.zarobki >= 0) & (df.zarobki<=95.45775)]
print('Pozostało ' + str(df.shape[0]) + ' obserwacji.')
4.1.2. Wzrost¶
Usuwam ze zbioru:
- zbyt wysokie kobiety,
- zbyt niskie kobiety,
- zbyt wysokich mężczyzn.
df = df[((df.plec=='mężczyzna') & (df.wzrost>=158)) | (((df.plec=='kobieta') & (df.wzrost>=148)) & ((df.plec=='kobieta') & (df.wzrost<=180)))]
print('Pozostało ' + str(df.shape[0]) + ' obserwacji.')
Po usunięciu wartości odstających, z początkowych 1379 obserwacji pozostało 1304. Wprawdzie próba została nieco uszczuplona, ale na szczęście nie na tyle by miało to negatywny wpływ na proces wnioskowania parametrycznego.
4.2. Normalizacja zmiennych numerycznych¶
Punkt zostaje pominięty. Przyczyna jest prosta: po znormalizowaniu danych numerycznych stracę możliwość interpretacji modelu regresji liniowej. Zatem nie będę w stanie powiedzieć ile warty jest każdy rok dodatkowej nauki ;) Pomijam ten krok.
5. Postawienie hipotez¶
Dane są przygotowane, można zatem przejść do ich głębszej analizy. Jej głownym celem jest postawienie kilku hipotez, które następnie będę badać. W procesie analizy i formułowania hipotez wykorzystam wizualizacje danych pochodzących z przygotowanej próbki.
Podczas formułowania hipotezy zerowej zakładał będę równość badanych parametrów, ew. brak zależności. Podobnie jak w przypadku spraw sądowych, gdzie na początku zakłada się niewinność oskarżonego i dopiero w dalszej kolejności zbiera materiał dowodowy który go obciąża. Jeśli liczba dowodów będzie odpowiednio duża, to odrzuca się hipotezę zerową i przyjmuje hipotezę alternatywną. Proste :)
5.1. Poziom zarobków a płeć¶
Pierwsza hipoteza będzie dotyczyć pytania, które przyszło mi ona do głowy, gdy pierwszy raz zobaczyłem dane: czy płeć danej osoby ma wpływ na zarobki?
summary = df.groupby('plec')['zarobki'].describe().sort_values(by='mean', ascending=False)
summary['median'] = df.groupby('plec')['zarobki'].median().astype(float)
summary
Pierwsze dowody ukazują, że poziom rocznych zarobków wśród kobiet i mężczyzn znacząco się różni. Z ostatecznym osądem jednak jeszcze się wstrzymam.
plt.figure(figsize=(15,7))
sns.set(font_scale=1.5)
sns.barplot('plec','zarobki', data=df, palette='Blues_d', capsize=0.03, order=summary.index).set(title='Zarobki per płeć')
plt.show()
Hipoteza 1: nie ma istotnej różnicy pomiędzy poziomem zarobków różnych płci.
5.2. Poziom zarobków wśród poszczególnych ras¶
summary = df.groupby('rasa')['zarobki'].describe().sort_values(by='mean', ascending=False)
summary['median'] = df.groupby('rasa')['zarobki'].median().astype(float)
summary
Dla zwiększenia czytelności posortowałem rasy po średniej wartości zarobków, w sposób malejący. Średnio najwięcej zarabiają osoby rasy białej. Co ciekawe, kiedy spojrzymy na medianę, to największe zarobki mają przedstawiciele (a raczej przedstawiciel) rasy czarnej. Odchylenie standardowe dla tej rasy jest również najmniejsze.
plt.figure(figsize=(15,7))
sns.set(font_scale=1.5)
sns.barplot('rasa','zarobki', data=df, palette='Blues_d', capsize=0.03, order=summary.index).set(title='Zarobki per rasa')
plt.show()
Powyższy wykres przedstawia średni poziom zarobkówo dla poszczególnych ras. Zgodnie z wynikami z tabeli podsumowującej, najwięcej zarabiają osoby rasy białej, a najmniej rasy hiszpańskiej. Czarne linie prezentują granice przedziałów ufności dla poszczególnych ras. W praktyce przedziały ufności informują o tym "na ile możemy ufać danej wartości".
Odnosząc ostatnie zdanie do wyników widocznych na wizualizacji można stwierdzić, że powinienem mieć najmniejsze zaufanie do średniego poziomu zarobków przedstawionych dla rasy "inna". "Szeroki" przedział mówi o tym, że rzeczywista wartość tego parametru (średnia zarobków dla rasy "inna") dla całej populacji (podczas gdy my posługujemy się jedynie próbą) mieści się w przedziale ok. [20k,36k]. Pewnie zauważyłeś, że im mniej obserwacji danej kategorii tym "szerszy" przedział. Wynika to m.in. z faktu, że im więcej obserwacji jest poddanych badaniu, tym większą część całej populacji jest zawarta w analizie i tym większa pewność do uzyskanych wyników.
Hipoteza 2: nie ma istotnej różnicy pomiędzy poziomem zarobów wśród różnych ras.
5.3. Poziom zarobków a rasa i płeć¶
Idę o krok dalej. Sprawdzę czy kobiety zawsze zarabiają mniej. Czy jakikolwiek wpływ na poziom rocznych zarobków ma rasa?
summary = df.groupby(['plec','rasa'])['zarobki'].describe().sort_values(by='mean', ascending=False)
summary['median'] = df.groupby(['plec','rasa'])['zarobki'].median().astype(float)
summary
Statystycznie najwięcej wśród mężczyzn zarabiają osoby rasy białej. Wśród kobiet najwięcej zarabiają osoby rasy 'inna'. To małe zaskoczenie. Należy jednak mieć na uwadze liczbę obserwacji należących do tej rasy: mamy zaledwie 9 mężczyzn i 17 kobiet. Może się to okazać nieco za mało by na tak małej próbie wnioskować o całej populacji.
plt.figure(figsize=(15,7))
sns.set(font_scale=1.5)
sns.barplot('plec', 'zarobki', hue='rasa', data=df, order = ['mężczyzna','kobieta'], palette='Blues_d', capsize=0.02).set(title='Zarobki per rasa i płeć', xlabel='płeć')
plt.show()
Z powyższego wykresu jasno wynika, że kobiety w każdej z ras zarabiają znacznie mniej od mężczyzn. Niestety nie ma żadnego przypadku w którym sytuacja jest odwrotna. Najbardziej zbliżony poziom zarobków jest w obserwacjach rasy 'inna'. Sprawdzę więc czy mam wystarczające dowody na to by stwierdzić, że mężczyźni tej rasy zarabiają istotnie więcej od kobiet.
Hipoteza 3: nie ma istotnej różnicy pomiędzy poziomem zarobków wśród kobiet i mężczyzn rasy 'inna'.
5.4. Poziom zarobków a wzrost¶
sns.set(font_scale=1.5)
sns.lmplot('wzrost', 'zarobki', data=df, palette='Blues', size=8, aspect=1.6).set(title='zarobki vs wzrost')
plt.show()
Hipoteza 4: nie istnieje zależność pomiędzy wzrostem a poziomem zarobków
5.5. Poziom zarobków a edukacja¶
W tym aspekcie chciałbym sprawdzić 2 rzeczy:
- Jak na poziom zarobków wpływał dodatkowy rok edukacji?
- Jak na poziom zarobków wpływał kolejny etap edukacji?
5.5.1. Poziom zarobków a lata edukacji¶
sns.set(font_scale=1.5)
sns.lmplot('lata_edukacji', 'zarobki', data=df, palette='Blues', size=8, aspect=1.6).set(title='lata_edukacji vs zarobki')
plt.show()
Unikalnych wartości w zmiennej 'lata_edukacji' jest relatywnie mało, więc pozwolę tu sobie na nieco głębszą analizę.
summary = df.groupby('lata_edukacji')['zarobki'].describe().sort_values(by='mean', ascending=False)
summary['median'] = df.groupby('lata_edukacji')['zarobki'].median().astype(float)
summary
No i mamy spore zaskoczenie. Wykres sugeruje dodatnią (ew. pozytywną) zależność liniową. Są jednak miejsca w których zależność jest odwrotna:
- 9 i 10 lat edukacji - 17.5k vs 15.1k rocznie,
- 14 i 15 lat edukacji - 31.6k vs 27.2k rocznie.
Na ten moment nie wiem z czego to może wynikać. Ciekaw jestem czy przy takich anomaliach lata edukacji wpływają w sposób istotny statystycznie na poziom zarobków.
Hipoteza 5: nie ma istotnej różnicy w poziomie zarobków pomiędzy osobami o różnej liczbie lat edukacji.
5.5.2. Poziom zarobków a zdobyte wykształcenie¶
summary = df.groupby('wyksztalcenie')['zarobki'].describe().sort_values(by='mean', ascending=False)
summary['median'] = df.groupby('wyksztalcenie')['zarobki'].median().astype(float)
summary
plt.figure(figsize=(15,7))
sns.set(font_scale=1.5)
sns.barplot('wyksztalcenie', 'zarobki', data=df, palette='Blues_d', capsize=0.02).set(title='Zarobki per poziom wykształcenia', xlabel='wykształcenie')
plt.show()
Hipoteza 6: nie ma istotnej różnicy w poziomie zarobków pomiędzy osobami o różnym poziomie wykształcenia.
5.6. Poziom zarobków a płeć i edukacja¶
summary = df.groupby(['plec','wyksztalcenie'])['zarobki'].describe()
summary['median'] = df.groupby(['plec','wyksztalcenie'])['zarobki'].median().astype(float)
summary
plt.figure(figsize=(15,7))
sns.set(font_scale=1.5)
sns.barplot('plec', 'zarobki', hue='wyksztalcenie', data=df, order = ['mężczyzna','kobieta'], palette='Blues_d', capsize=0.02).set(title='Zarobki per poziom wykształcenia i płeć', xlabel='płeć')
plt.show()
Spora niespodzianka. Wszystko wskazuje na to, że kobiety o wykształceniu wyższym zarabiają mniej niż mężczyźni o wykształceniu średnim. Tu chyba mamy sprawę rozstrzygniętą. Co zatem ma większy wpływ na poziom zarobków: płeć, czy wykształcenie?
Hipoteza 7: wykształcenie nie ma większego wpływu na poziom zarobków niż płeć.
6. Badanie zależności pomiędzy zmiennymi¶
W tym punkcie będę badać korelacje i zależność pomiędzy poszczególnymi zmiennymi. Metod badania jest kilka i można je podzielić na 3 podstawowe grupy:
- liczbowa zależność/korelacja (określająca siłę powiązania pomiędzy badanymi zmiennymi),
- test zależności pomiędzy zmiennymi (odpowiada na pytanie czy istnieje statystyczna zależność pomiędzy badanymi zmiennymi)
- interpretacja wizualizacji.
Pytanie jednak, po co w ogóle badać zależność pomiędzy zmiennymi? Otóż, zakładałem od początku, że do testowania hipotez zastosuję wnioskowanie parametryczne oparte o model regresji liniowej. Model ten posiada jedno dosyć poważne założenie: brak współliniowości pomiędzy zmiennymi objaśniającymi. Muszę zweryfikować, czy nie mam w zbiorze zmiennych, które potencjalnie mogą być przyczyną współliniowości.
6.1. Badanie korelacji (zmienne numeryczne)¶
Do badania korelacji pomiędzy zmiennymi numerycznymi można użyć dówch podstawowych metod:
- Korelacji Pearsona.
- Korelacji rang Spearmana.
Pierwsza z nich ma dwa ograniczenia:
- jej użycie jest sensowne jedynie dla zmiennych o rozkładzie normalnym,
- jest wrażliwa na wartości odstające.
Żadna ze zmiennych numerycznych, które są widoczne w próbie nie cechuje się rozkładem normalnym. Korelacja Pearsona odpada. Na szczęcie oba powyższe ograniczenia można "obejść" dzięku zastosowaniu korelacji rang Spearmana. Obserwacje odstające stają się w niej nieistotne, gdyż zamiast wartości zmiennoprzecinkowych kolejne wartości są zamieniane na rangi.
Do badania korelacji zmiennych numerycznych włączę płeć, gdyż jako zmienna kategoryczna przyjmuje ona jedynie dwie wartości (mężczyzna, kobieta). Znając założenia korelacji ran Spermana, wiem że mogę sobie na to pozwolić :)
Przedstawienie wyników:¶
col = ['zarobki','wzrost', 'wiek', 'lata_edukacji', 'plec']
corr = pd.DataFrame(stats.spearmanr(df[col])[0], columns=col, index=col)
corr
Wnioski:
- Na pierwszy rzut oka widać, że mogę mieć problem z ujęciem w jednym modelu zmiennych: 'plec', 'wzrost'. Są one ze sobą mocno skorelowane i umieszczenie ich w jednym modelu nie jest dobrym pomysłem.
- Najmocniej skorelowaną zmienną ze zmienną celu ('zarobki') jest 'plec'.
- Na drugim miejscu jest 'wzrost' (skoro wiem, że najmocniej skorelowana jest 'plec', która jest bardzo mocno skorelowana ze zmienną 'wzrost', to raczej nie będę jej umieszczać w moim modelu. Jesli wystąpi współliniowość, to będzie ona pierwszym podejrzanym).
- Na trzecim 'lata_edukacji' (tu niewielkie zaskoczenie bo istnieje mocne przekonanie, że to właśnie od czasu poświęconego na edukację zależą nasze przyszłe zawodowe sukcesy, mierzone poziomem zarobków).
6.2. Badanie zależności (zmienne numeryczne vs zmienne kategoryczne)¶
W przypadku zmiennych numerycznych i kategorycznych nie mówimy już o korelacji, lecz o poziomie zależności. Miarą, która jest najbliższa współczynnikowi korelacji, w tym przypadku jest współczynnik korelacji wielorakiej (współczynnik R w modelach regresji) – użyję jej do badania zależności pomiędzy zmiennymi kategorycznymi (nominalnymi i porządkowymi) i numerycznymi (liczby całkowite i zmiennoprzecinkowe).
df.dtypes
Zarobki vs (płeć, rasa, wykształcenie):¶
# płeć
model = sm.ols(formula='zarobki~plec',data=df)
zar_ple = np.sqrt(model.fit().rsquared)
# rasa
model = sm.ols(formula='zarobki~rasa',data=df)
zar_ras = np.sqrt(model.fit().rsquared)
# wykształcenie
model = sm.ols(formula='zarobki~wyksztalcenie',data=df)
zar_wyk = np.sqrt(model.fit().rsquared)
zar = pd.Series([zar_ple, zar_ras, zar_wyk])
Wzrost vs (płeć, rasa, wykształcenie):¶
# płeć
model = sm.ols(formula='wzrost~plec',data=df)
wzr_ple = np.sqrt(model.fit().rsquared)
# rasa
model = sm.ols(formula='wzrost~rasa',data=df)
wzr_ras = np.sqrt(model.fit().rsquared)
# wykształcenie
model = sm.ols(formula='wzrost~wyksztalcenie',data=df)
wzr_wyk = np.sqrt(model.fit().rsquared)
wzr = pd.Series([wzr_ple, wzr_ras, wzr_wyk])
Lata edukacji vs (płeć, rasa, wykształcenie):¶
# płeć
model = sm.ols(formula='lata_edukacji~plec',data=df)
lat_ple = np.sqrt(model.fit().rsquared)
# rasa
model = sm.ols(formula='lata_edukacji~rasa',data=df)
lat_ras = np.sqrt(model.fit().rsquared)
# wykształcenie
model = sm.ols(formula='lata_edukacji~wyksztalcenie',data=df)
lat_wyk = np.sqrt(model.fit().rsquared)
lat = pd.Series([lat_ple, lat_ras, lat_wyk])
Wiek vs (płeć, rasa, wykształcenie):¶
# płeć
model = sm.ols(formula='wiek~plec',data=df)
wie_ple = np.sqrt(model.fit().rsquared)
# rasa
model = sm.ols(formula='wiek~rasa',data=df)
wie_ras = np.sqrt(model.fit().rsquared)
# wykształcenie
model = sm.ols(formula='wiek~wyksztalcenie',data=df)
wie_wyk = np.sqrt(model.fit().rsquared)
wie = pd.Series([wie_ple, wie_ras, wie_wyk])
Podsumowanie:¶
summary = pd.DataFrame([zar, wzr, lat, wie])
summary.index = ['zarobki','wzrost','lata_edukacji','wiek']
summary.columns = ['plec','rasa','wyksztalcenie']
summary
Wnioski:
Największy poziom zależności jest widoczny pomiędzy zmiennymi: 'lata_edukacji', 'wyksztalcenie'. Nie jest to dziwne, gdyż zmienna 'wyksztalcenie' została utworzona na podstawie zmiennej 'lata_edukacji' i dodana do zbioru przeze mnie.
Druga w kolejności jest zalezność pomiędzy zmiennymi: 'wzrost', 'plec'. O tym również pisałem wcześniej i chyba nie trzeba dowodzić tej zależności:)
Trzecią parą o najwyższym poziomie zalezności są: 'zarobki', 'plec'. Mam więc kolejny dowód na to, że zarobki są zależne od płci.
Dopiero na czwartym miejscu jest para: 'zarobki', 'wyksztalcenie'. Wydawać by się mogło, że większy wpływ na poziom zarobków ma wykształcenie niż płeć. To oczywiście nie jest jeszcze ostateczny dowód, lecz kolejna przesłanka.
Pozostałe zależności nie są aż tak istotne.
6.3. Badanie zależności (zmienne kategoryczne)¶
Do badania zależności pomiędzy poszczególnymi zmiennymi kategorycznymi użyję współczynnika V Craméra. Nie chcę go szczegółowo omawiać, ani opisuwać czemu akurat zdecydowałem się na użycie tej miary. Być może pochylę się nad nim w jendym z kolejnych wpisów. W tym momencie warto zaznaczyć, że większość znanych mi metod badania zależności pomiędzy zmmiennymi kategorycznymi bazuje na wynikach tabeli krzyżowej, przedstawiającej liczbę obserwacji spełniających wybrane kryteria.
Niestety nie udało mi się znaleźć żadnej biblioteki z gotową metodą wyznaczania współczynnika V Craméra, dlatego muszę sam ją zaimplementować.
def CramersV(tab):
a = stats.chi2_contingency(tab)[0]/sum(tab.sum())
b = min(tab.shape[0]-1, tab.shape[1]-1,)
return(np.sqrt(a/b))
def CalculateCrammersV(tab):
ret = []
for m in tab:
row = []
for n in tab:
cross_tab = pd.crosstab(tab[m],tab[n])
row.append(CramersV(cross_tab))
ret.append(row)
return pd.DataFrame(ret, columns=tab.columns, index=tab.columns)
Podsumowanie:¶
CalculateCrammersV(df[['plec','rasa','wyksztalcenie']])
Wnioski:
Najwyższy poziom zależności występuje pomiędzy zmiennymi: 'rasa', 'wyksztalcenie'. Nie mniej, jest to ciągle bardzo mała zależność, którą w zasadzie można pominąć.
7. Weryfikacja hipotez¶
Podejście do weryfikacji hipotez¶
- Na potrzeby weryfikacji hipotez zakładam poziom istotności alfa=0.05.
- Dla każdej hipotezy wyznaczam hipotezę zerową (H0) i hipotezę alternatywną (H1).
Do weryfikowania każdej z nich, mógłbym podejść z osobna, budując oddzielne modele. To jendak zwiększyłoby szanse zaobserwowania zjawiska nazywanego paradoksem Simpsona.
Paradoks Simpsona¶
W statystyce jest to zjawisko opisujące zmianę wyników wnioskowania po uwzględnieniu nowych informacji/zmiennych. Dla przykładu, budując model w którym użyję jedynie dwóch zmiennych: 'zarobki' i 'wzrost', nie wiem czy na tą relację nie mają wpływu inne zmienne. Dodanie kolejnych zmiennych do analizy może (ale nie musi) wpłynąć na wynik analizy.
By dokładnie pokazać jak wygląda powyższe zjawisko, zbuduję dwa modele (Model 1 i Model 2) do weryfikacji hipotezy 4. Mówi ona o tym, że nie istnieje zależność pomiędzy wzrostem a poziomem zarobków.
Model 1¶
# Buduję model regresji liniowej
model = sm.ols(formula='zarobki~wzrost', data=df)
# Drukuję podsumowanie modelu
print(model.fit().summary())
Interpretacja wyników:¶
- Wartość Prob (F-statistic) jest mniejsza od zakładanego poziomu istotności, dlatego przyjmuję, że model jako całość jest istotny statystycznie.
- Komunikat [2] sugeruje, że w modelu zachodzi zjawisko współliniowości, lub inny problem numeryczny. Na potrzeby tego eksperymentu ignoruję ten komunikat. O problemie ze zmienną 'wzrost' wspominałem w poprzednich punktach.
- Wartość 'coef' dla Intercept, to wartość współczynnika beta_0 (wyraz wolny) z modelu regresji liniowej.
- Wartość 'coef' dla 'wzrost', to wartość współczynnika beta_1 z modelu regresji liniowej.
- Wartości P>|t|, to p-value dla poszczególnych zmiennych. Zarówno dla 'wzrostu' jak i dla wyrazu wolnego są mniejsze niż alfa. W modelu regresji liniowej należy je interpretować w następujący sposób: p-value oznacza prawdopodobieństwo, że rzeczywista wartość współczynnika beta (dla całej populacji, a nie tylko dla próby) wynosi zero, przy zaobserwowanym poziomie współczynnika beta. Czemu porównuję do zera? Ponieważ w modelu regresji liniowej, każda zmienna której współczynnik wynosi zero (ew. istnieje wysokie prawdopodobieństwo, że może wynosić zero) jest zbędna, gdyż nic nie wnosi. W żaden sposób nie przyczynia się ona do lepszego dopasowania modelu do zmiennej celu. Chcemy unikać takich sytuacji, więc interesować mnie będą jedynie te zmienne, dla których wartość p-value jest możliwie jak najmniejsza (prawdopodobieństwo popełnienia błędu I rodzaju jest niewielkie).
- Wartości p-value dla poszczególnych współczynników są bardzo małe, a więc można podejrzewać, że wzrost istotnie wpływa na zarobki.
- Wartości: [0.025 0.975] - to przedziały ufności na poziomie 95%, dla danej zmiennej. Oznaczają one, że z przyjętym poziomem ufności równym 95%, prawdziwa wartośc współczynnika leży w wyznaczonym przedziale, np. zgodnie z powyższym modelem i przyjętym poziomem ufności 95%, prawdziwa wartość współczynnika beta dla całej populacji, dla zmiennej 'wzrost' leży w przedziale [0.682, 0.917].
- Zgodnie z wynikami modelu, estymacja rocznych zarobków dla osoby, która ma 175 cm wzrostu wynosi: $\$$33.71k. Wynika to z prostego równania: zarobki = Intercept + wzrost * 0.7994.
- Każdy centymetr wzrostu jest wart ok. $\$$799 do rocznej pensji. Ale czy na pewno? By to sprawdzić odaję do modelu kolejną zmienną: 'plec'.
Model 2¶
model = sm.ols(formula='zarobki~wzrost+plec', data=df)
print(model.fit().summary())
Interpretacja wyników:¶
- Komunikat [2] informuje o problemie ze współliniowością. Jeszcze raz ignoruję go na potrzeby eksperymentu :)
- Model jako całość jest istotny statystycznie (Prob (F-statistic): < alfa).
- Wartości parametrów dla zmiennych 'wzrost' i 'plec' również są istotne statystycznie.
- Najmniejszą pewność mam co do współczynnika beta_0 (p-value dla Intercept wynosi 0.011. Prawdopodobieństwo popełnienia błędu w tym przypadku, to ok 1.1%).
- Zgodnie z wynikami tego modelu, mężczyźni zarabiają ok. $\$$12.18k. rocznie więcej niż kobiety.
- Każdy centymetr wzrostu jest już wart wart tylko ok. $\$$355.
Po dodaniu tylko jednej zmiennej do pierwszego modelu, 1 cm wzrostu jest już wart dwukrotnie mniej. To jest właśnie paradoks Simpsona. Im więcej zmiennych o wysokiej jakości (takich, które nie wnoszą do modelu współliniowości, o możliwie niskim p-value i współczynniku beta różnym od zera) w modelu, tym mniejsze prawdopodobieństwo wystąpienia paradoksu Simpsona.
Kilka dodatkowych uwag:¶
- Wcześniej podjęte kroki (m.in. badanie korelacji) pokazały, że zmienna 'wzrost' powinna zostać pominięta. Powoduje ona powstawanie współliniowości, dlatego pominę ją w dalszych rozważaniach. Jednocześnie nie jestem w stanie powiedzieć nic o wpływie wzrostu na poziom zarobków. Oczywiście, istnieje wysoka korelacja pomiędzy zmiennymi, natomiast zagłębiając się w dane można stwierdzić, że w rzeczywistości wynika ona z czego innego (wpływ płci na zarobki przy jednoczesnej wysokiej zależności pomiędzy płcią a wzrostem). Hipotezę 4 (dotyczącą wpływu wzrostu na zarobki) wykreślam zatem z dalszych rozważań, gdyż niezależnie od uzyskanych wyników, nie będą one istotne statystycznie.
- Nie powinienem dodawać do jednego modelu zmiennych: 'wyksztalcenie', 'lata_edukacji'. Są one od siebie mocno zależne.
- We wnioskowaniu parametrycznym mówi się o całej populacji badając jedynie jej wycinek. Zatem nawet, gdy wartości średnie sugerowały różnicę w poziomie zarobków osób różnych ras, to wyniki wnioskowania parametrycznego wcale nie muszą tego potwierdzać.
Model 3¶
model = sm.ols(formula='zarobki~plec+rasa+lata_edukacji+wiek', data=df)
print(model.fit().summary())
Interpretacja wyników:¶
Model jako całość jest istotny statystycznie (Prob (F-statistic) < alfa). Mogę zatem przejść do weryfikacji hipotez w oparciu o ten model.
Weryfikacja hipotez:¶
Hipoteza 1.
H0 - nie ma istotnej różnicy pomiędzy poziomem zarobków różnych płci.
H1 - jest istotna różnica pomiędzy poziomem zarobków różnych płci.
Spośród dwóch osób o różnej płci, a tym samym wzroście, mężczyzna zarabia o $\$$17.01k rocznie więcej niż kobieta. Wynik jest istotny statystycznie. Z maksymalnym prawdopodobieństwem popełnienia błędu równym 5%, odrzucam H0 i przyjmuję H1: jest istotna różnica pomiędzy poziomem zarobków różnych płci.
Hipoteza 2.
H0 - nie ma istotnej różnicy pomiędzy poziomem zarobów wśród różnych ras.
H1 - jest istotna różnica pomiędzy poziomem zarobków wśród różnych ras.
Jako, że p-value dla rasy: czarnej, hiszpańskiej i innej jest większe niż alfa, to nie pozwala mi ti przypuszczać, że prawdziwa wartość współczynnika beta dla nich jest różna od zera. Nie mogę zatem potwierdzić, że istnieje istotna statystycznie różnica pomiędzy śrendim poziomem zarobków osób różnych ras. Przyjmuję H0 - nie ma istotnej różnicy pomiędzy poziomem zarobków wśród różnych ras.
Hipoteza 5.
H0 - nie ma istotnej różnicy w poziomie zarobków pomiędzy osobami o różnej liczbie lat edukacji.
H1 - jest istotna różnica w poziomie zarobków pomiędzy osobami o różnej liczbie lat edukacji.
Niezależnie od płci i wieku, każdy dodatkowy rok edukacji sprawia, że dana osoba zarabia średnio o $\$$3.07k rocznie więcej. Wynik jest istotny statystycznie Zatem warto się uczyć :) Z maksymalnym prawdopodobieństwem popełnienia błędu równym 5%, odrzucam H0 i przyjmuję H1: jest istotna różnica w poziomie zarobków pomiędzy osobami o różnej liczbie lat edukacji.
Model 4¶
W celu weryfikacji hipotezy 3, muszę zbudować osobny model i zbadać tzw. efekt oddziaływania. Pierwsze badania wskazywały na to, że w przypadku rasy 'inna' efekt oddziaływania płci na zarobki może być inny niż w przypadku pozostałych ras. Przesłankami ku temu były małe różnice w zarobkach obu płci dla tej rasy.
model = sm.ols(formula='zarobki~plec:rasa+wiek+lata_edukacji', data=df)
print(model.fit().summary())
Interpretacja wyników:¶
Model jako całość jest istotny statystycznie (Prob (F-statistic) < alfa). Mogę zatem przejść do weryfikacji hipotez w oparciu o ten model.
Weryfikacja hipotezy 3:¶
H0 - nie ma istotnej różnicy pomiędzy poziomem zarobków wśród kobiet i mężczyzn rasy 'inna',
H1 - jest istotna różnica pomiędzy poziomem zarobków wśród kobiet i mężczyzn rasy 'inna'.
Zgodnie z wynikami modelu, nie można stwierdzić czy istnieje statystycznie istotna różnica pomiędzy poziomem zarobów wśród kobiet i mężczyzn rasy 'inna'. Wartość p-value w obu przypadkach jest większa od zakładanego poziomu istotności alfa. Przyjmuję H0 - nie ma istotnej różnicy pomiędzy poziomem zarobków wśród kobiet i mężczyzn rasy 'inna'.
Model 5¶
model = sm.ols(formula='zarobki~plec+rasa+wiek+wyksztalcenie', data=df)
print(model.fit().summary())
Interpretacja wyników:¶
- Model jako całość jest istotny statystycznie (Prob (F-statistic) < alfa). Mogę zatem przejść do weryfikacji hipotez w oparciu o ten model.
- Każdy kolejny etap edukacji niesie za sobą wyższe roczne zarobki.
- Największy przeskok w zarobkach jest pomiędzy poziomem szkołu podstawowej, a szkoły średniej.
Weryfikacja hipotezy 6:¶
H0 - nie ma istotnej różnicy w poziomie zarobków pomiędzy osobami o różnym poziomie wykształcenia.
H1 - jest istotna różnica w poziomie zarobków pomiędzy osobami o różnym poziomie wykształcenia.
Zgodnie z wynikami modelu, można stwierdzić, że istnieje statystycznie istotna różnica pomiędzy poziomem zarobów wśród osób o różnym poziomie wykształcenia. We wszystkich przypadkach wartość p-value jest mniejsza od zakładanego poziomu istotności alfa. Odrzucam H0 i przyjmuję H1 - jest istotna różnica w poziomie zarobków pomiędzy osobami o różnym poziomie wykształcenia.
Weryfikacja hipotezy 7:¶
H0 - wykształcenie nie ma większego wpływ na poziom zarobków niż płeć.
H1 - wykształcenie ma większy wpływ na poziom zarobków niż płeć.
Kończąc studia dana osoba, niezależnie od rasy, płci i wieku zarabia o $\$$26.663k rocznie więcej niż osoba która ma wykształcenie podstawowe. Jesli chodzi o płeć, to różnica w zarobkach u kobiet i mężczyzn tej samej rasy, wieku i wykształcenia wynosi $/$$17.189k rocznie. We wszystkich przypadkach wartość p-value jest mniejsza od zakładanego poziomu istotności alfa. Odrzucam H0 i przyjmuję H1 - wykształcenie ma większy wpływ na poziom zarobków niż płeć.
8. Podsumowanie¶
Wyniki projektu w skrócie:
- płeć ma wpływ na zarobki,
- wykształcenie ma większy wpływ na zarobki niż płeć,
- nie ma istotnej różnicy pomiędzy poziomem zarobków wśród różnych ras,
- wzrost nie ma wpływu na poziom zarobków.