Kernel

1 Analiza opisowa

Jednym ze sposobów zrozumienia, jak działa rząd miasta, jest spojrzenie na to, kogo zatrudnia i jak jego pracownicy są wynagradzani. Dane te zawierają nazwiska, nazwę stanowiska i wynagrodzenie pracowników miasta San Francisco w ujęciu rocznym od 2011 do 2014 roku.

Oto kilka pomysłów na eksplorację danych:

• Jak zmieniały się wynagrodzenia w czasie między różnymi grupami ludzi?

• Jak płaca podstawowa, wynagrodzenie za nadgodziny i świadczenia są rozdzielane pomiędzy różne grupy?

• Czy w tym zestawie danych istnieją dowody na dyskryminację płacową ze względu na płeć?

• Jak przydzielany jest budżet w zależności od grupy i zakresu obowiązków?

# wymiary ramki:
dim(salaries)

## [1] 148654     13

# nazwy kolumn:
names(salaries)

##  [1] "Id"               "EmployeeName"     "JobTitle"         "BasePay"         
##  [5] "OvertimePay"      "OtherPay"         "Benefits"         "TotalPay"        
##  [9] "TotalPayBenefits" "Year"             "Notes"            "Agency"          
## [13] "Status"

1.1 Histogramy

hist(salaries$TotalPay,main="Total Pay", xlab="Pay (in dollars)")
abline(v = mean(salaries$TotalPay),lty="dashed")
abline(v = median(salaries$TotalPay))
legend("topright", legend=c("Mediana","Średnia"),lty=c("solid","dashed"))

par(mfrow=c(2,2))
hist(salaries$TotalPay,main="Total Pay, default breaks", xlab="Pay (in dollars)")
hist(salaries$TotalPay,main="Total Pay, breaks=100", xlab="Pay (in dollars)", breaks=100)
hist(salaries$TotalPay,main="Total Pay, breaks=1000", xlab="Pay (in dollars)",breaks=1000)

hist(salaries$TotalPay,main="Total Pay, Zoomed-in", xlab="Pay (in dollars)", xlim=c(0,1e5), breaks=1000)

salaries2 <- subset(salaries, JobTitle=="Firefighter" & Status=="FT")
dim(salaries2)

## [1] 738  13

par(mfrow=c(2,2))
hist(salaries2$TotalPay,main="Firefighters, default breaks", xlab="Pay (in dollars)")
hist(salaries2$TotalPay,main="Firefighters, breaks=30", xlab="Pay (in dollars)", breaks=30)
hist(salaries2$TotalPay,main="Firefighters, breaks=100", xlab="Pay (in dollars)", breaks=100)
hist(salaries2$TotalPay,main="Firefighters, breaks=1000", xlab="Pay (in dollars)",breaks=1000)

1.2 Wykresy pudełkowe

par(mfrow=c(1,1))
boxplot(salaries$TotalPay,main="Total Pay, breaks=1000", ylab="Pay (in dollars)")

2 Estymacja funkcji gęstości

Pierwszy raport dotyczy nieparametrycznej estymacji gęstości. Klasycznym nieparametrycznym estymatorem gęstości jest histogram, który dostarcza nieciągłe i stałe oszacowania. W tym raporcie skupiono się na niektórych alternatywach, które zapewniają ciągłe lub nawet gładkie oszacowania zamiast.

Metody kernelowe stanowią ważną klasę gładkich estymatorów gęstości i zaimplementowane są przez funkcję R density(). Estymatory te są w zasadzie tylko lokalnie ważonymi średnimi, a ich obliczenie jest stosunkowo proste w teorii. W praktyce, różne wybory sposobu implementacji obliczeń mogą jednak mieć duży wpływ na rzeczywisty czas obliczeń, a implementację kernelowych estymatorów gęstości zilustruje trzy punkty:

• jeśli to możliwe, wybierz wektoryzowane implementacje w R,
• jeśli niewielka strata w dokładności jest do zaakceptowania, przybliżone rozwiązanie może być o rzędy wielkości szybsze niż implementacja literalna,
• czas potrzebny do numerycznej oceny różnych funkcje elementarne może bardzo zależeć od funkcji i sposobu implementacji obliczeń.

Metody kernelowe opierają się na jednym lub więcej parametrach regularności, które muszą być dobrane tak, aby osiągnąć właściwą równowagę w dostosowaniu do danych bez zbytniego dostosowywania się do losowej zmienności w danych.

Wybór odpowiedniej ilości regularności jest równie ważny jak wybór metody do użycia w pierwszej kolejności. W rzeczywistości może być ważniejszy. Tak naprawdę nie mamy kompletnej implementacji nieparametrycznego estymatora dopóki nie zaimplementujemy automatycznego, opartego na danych sposobu wyboru ilości regulacji.

Implementacja tylko obliczeń dla oceny estymatora jądra, powiedzmy, i pozostawiając to całkowicie użytkownikowi wyboru szerokości pasma jest pracą w połowie wykonaną. Metody i implementacje do wyboru szerokości pasma są więc w tym raporcie omówione dość szczegółowo.

W ostatniej części przeprowadzona jest analiza prawdopodobieństwa. Robi się to w celu dalszego wyjaśnienia, dlaczego potrzebne są estymatory z regularyzacją w celu uniknięcia nadmiernego dopasowania do danych, oraz dlaczego nie istnieje w ogóle nieparametryczny maksymalnego prawdopodobieństwa estymatora gęstości. Regularyzację prawdopodobieństwamożna osiągnąć poprzez ograniczenie szacunków gęstości do rodziny coraz bardziej elastycznych gęstości parametrycznych, które są dopasowane do danych. Jest to znane jako metoda sit. Inne podejście opiera się na rozszerzeniach bazowych, ale w obu przypadkach automatyczny wybór wielkości regularności jest tak samo ważny jak w przypadku metod jądrowych.

Aby utworzyć wykres gęstości jądra, musisz oszacować gęstość jądra. W tym celu można użyć funkcji density, a następnie przekazać obiekt density do funkcji plot.

# dane
set.seed(14012021)
data <- rnorm(200, mean = 4)

# Kernel density estimation
d <- density(data)

# Kernel density plot
plot(d, lwd = 2, main = "Default kernel density plot")

Argument jądra funkcji gęstości domyślnie używa jądra gaussowskiego (kernel = “gaussian”), ale dostępnych jest więcej typów jądra, takich jak “prostokątne”, “trójkątne”, “epanechnikov”, “biweight”, “cosine” i “optcosine”. Wybór będzie zależał od twoich danych, ale w większości scenariuszy wartość domyślna jest najbardziej zalecana.

# Data
set.seed(14012021)
data <- rnorm(200, mean = 4)

# Kernel density estimation
d <- density(data,
             kernel = "rectangular")

# Kernel density plot
plot(d, lwd = 2, main = "Rectangular kernel")

# Data
set.seed(14012021)
data <- rnorm(200, mean = 4)

# Kernel density estimation
d <- density(data,
             kernel = "triangular")

# Kernel density plot
plot(d, lwd = 2, main = "Triangular kernel")

# Data
set.seed(14012021)
data <- rnorm(200, mean = 4)

# Kernel density estimation
d <- density(data,
             kernel = "epanechnikov")

# Kernel density plot
plot(d, lwd = 2, main = "Epanechnikov kernel")

# Data
set.seed(14012021)
data <- rnorm(200, mean = 4)

# Kernel density estimation
d <- density(data,
             kernel = "biweight")

# Kernel density plot
plot(d, lwd = 2, main = "Biweight kernel")

# Data
set.seed(14012021)
data <- rnorm(200, mean = 4)

# Kernel density estimation
d <- density(data,
             kernel = "cosine")

# Kernel density plot
plot(d, lwd = 2, main = "Cosine kernel")

2.1 Selekcja pasma

Argument bw funkcji gęstości pozwala na zmianę używanego pasma wygładzania. Możesz przekazać wartość lub ciąg znaków podający regułę wyboru lub funkcję. Domyślną wartością jest “nrd0” (lub bw.nrd0(.)), która implementuje podejście oparte na zasadzie reguły kciuka :-) Inne dostępne opcje to:

2.1.1 Reguła Scotta (1992)

# Data
set.seed(14012021)
data <- rnorm(200, mean = 4)

# Kernel density estimation
d <- density(data,
             bw = "nrd")

# Kernel density plot
plot(d, lwd = 2, main = "nrd bandwidth")

2.1.2 Nieobciążona cross-walidacja

# Data
set.seed(14012021)
data <- rnorm(200, mean = 4)

# Kernel density estimation
d <- density(data,
             bw = "ucv")

# Kernel density plot
plot(d, lwd = 2, main = "ucv bandwidth")

2.1.3 Obciążona cross-walidacja

# Data
set.seed(14012021)
data <- rnorm(200, mean = 4)

# Kernel density estimation
d <- density(data,
             bw = "bcv")

# Kernel density plot
plot(d, lwd = 2, main = "bcv bandwidth") 

2.1.4 Metoda Sheather & Jones (1991)

# Data
set.seed(14012021)
data <- rnorm(200, mean = 4)

# Kernel density estimation
d <- density(data,
             bw = "SJ")

# Kernel density plot
plot(d, lwd = 2, main = "SJ bandwidth")

Ostrzeżenie!

Szerokość pasma musi być bardzo starannie dobrana! Mała szerokość pasma spowoduje powstanie nadmiernie dopasowanej krzywej, natomiast zbyt duża szerokość pasma spowoduje powstanie krzywej nadmiernie wygładzonej.

3 Ćwiczenie 1.

Uruchom demo estymatora funkcji gęstości kernel. Zmieniaj zarówno dane wejściowe, jak i opcje estymatora - szerokość pasma oraz rodzaj funkcji jądrowej. Czy widzisz istotne różnice w oszacowaniu?

#install.packages("remotes") #tylko raz! potem #
#install_github("hericks/KDE") #tylko raz! potem #
#install.packages("htmltools")
#remotes::install_github("hericks/KDE")
#install.packages("kde")
library(ggplot2)
library(KDE)
shiny_kde()

Shiny applications not supported in static R Markdown documents

Rzeczywiście można zauważyć wyraźne różnice. Zwiększenie wartości odchylenia standardowego powoduje, że rozkład staje się bardziej płaski i rozciągnięty, co oznacza większe rozproszenie danych wokół średniej. Z kolei zmniejszenie odchylenia standardowego sprawia, że rozkład staje się węższy i bardziej skupiony. Zmiana wartości średniej powoduje przesunięcie środka rozkładu na osi poziomej (oś x). Innymi słowy, cały rozkład przesuwa się w prawo lub w lewo w zależności od tego, czy średnia została zwiększona, czy zmniejszona, ale sam kształt rozkładu pozostaje taki sam. Natomiast modyfikacja szerokości pasma (bandwidth) wpływa na stopień, w jakim krzywa dopasowuje się do danych. Im mniejsza wartość pasma, tym bardziej krzywa staje się szczegółowa i lepiej odwzorowuje lokalne zmiany w danych. Jednak większa szerokość pasma powoduje wygładzenie krzywej, co prowadzi do bardziej ogólnego dopasowania i pominięcia drobnych szczegółów.

Jeśli chodzi o różnice pomiędzy argumentami jądra funkcji gęstości, to możemy zauwazyć, że rodzaj ‘rectangular’ jest najbardziej dopasowany do danych, można by stwierdzić, że aż zbyt dobrze dopasowany, ponieważ praktycznie wogóle nie wygładza rozkładu. Rodzaj ‘Logistic’, jest najbardziej “spłaszczony” i wygładzony ze wszystkich rozkładów. ‘Triweight’,‘Triangular’,‘Tricube’, dobrze dopasowują się do danych, ale niezbyt wygładzają rozkład. Posobnie ‘Epanechnikov’,‘Biweight’ i ‘Cosine’. ‘Silverman’ i ‘Sigmoid’ mają cechy zblizone do ‘Logistic’, lecz są one mniej spłaszczone. Choć jeszcze mnie spłaszczony jest klasyczny rozkład ‘Gaussian’.

4 Ćwiczenie 2.

Wykorzystując dowolną funkcję R do estymacji funkcji gęstości oszacuj jej przebieg dla wynagrodzeń (zbiór danych salaries) strażaków w San Francisco. Wykorzystaj metody graficzne dostępne w pakiecie ggplot2. Mile widziane przekroje oraz odpowiedzi na pytania badawcze zadane na wstępie.

salariesFFSanF <- subset(salaries, JobTitle=="Firefighter" & Agency=="San Francisco")
dim(salariesFFSanF)

## [1] 2359   13

4.1 Zaczniemy od pytań badawczych.

4.2 Jak zmieniały się wynagrodzenia w czasie między różnymi grupami ludzi?

4.3 Histogram wynagrodzeń dla strażaków z San Francisco.

sr_mediana_linie <- data.frame(nazwa = c("Średnia",
                                     "Mediana"),
                       wartość     = c(mean(salariesFFSanF$TotalPay),
                                  median(salariesFFSanF$TotalPay)))

ggplot(salariesFFSanF,aes(x = TotalPay)) +
  geom_histogram(binwidth=10000, color="purple",fill="yellow") +
  labs(x = "Pay (in dollars)", y = "Frequency", title = "Histogram wynagrodzeń strażaków w San Francisco") +
  geom_vline(data=sr_mediana_linie,aes(xintercept = wartość,
                               linetype = nazwa,
                               col = nazwa),linewidth=1,color="blue")+
  theme_minimal()

Możemy zauważyć, że badany rozkład wynagrodzeń jest prawoskośny i więcej strażaków zarabia ponad średnie wynagrodzenie dla ich stanowiska. Mamy również drugą, małą, lokalną modę w okolicach 0:50 000 USD, która zaniża nam wartość mediany.

4.4 Poniżej histogram zmian wynagrodzeń strażaków z SF na przestrzeni lat

sr_mediana_linie_2012 <- data.frame(nazwa = c("Średnia z 2012",
                                     "Mediana z 2012"), wartość = c(mean(salariesFFSanF$TotalPay[salariesFFSanF$Year==2012]),
 median(salariesFFSanF$TotalPay[salariesFFSanF$Year==2012])))

ggplot(salariesFFSanF,aes(x = TotalPay)) +
  geom_histogram(binwidth=25000, color="purple",fill="yellow") +
  labs(x = "Pay (in dollars)", y = "Frequency", title = "Histogram wynagrodzeń strażaków w San Francisco") +
  geom_vline(data=sr_mediana_linie,aes(xintercept = wartość,
                               linetype = nazwa,
                               col = nazwa),linewidth=1,color="blue")+
  theme_minimal()+
  facet_grid(. ~ Year)

Możemy zauważyć, że w porównaniu do 2012 roku w 2013 roku liczba strażaków, którzy zarabiali kwotę odpowiadającą średniej i medianie z 2012 roku, wzrosła. Jednak w 2014 roku wartości te powróciły do poziomów zbliżonych do tych z 2012 roku. Świadczy o tym zmiana średniego wynagrodzenia w analizowanych latach – po wzroście w 2013 roku osiągnęło ono najwyższy poziom, by następnie w 2014 roku spaść poniżej wartości odnotowanej w 2012 roku.

sr_mediany_lata <- data.frame(Rok = c(2012,2013,2014), 
                                    średnie = c(mean(salariesFFSanF$TotalPay[salariesFFSanF$Year==2012]),
mean(salariesFFSanF$TotalPay[salariesFFSanF$Year==2013]),
mean(salariesFFSanF$TotalPay[salariesFFSanF$Year==2014])),
                                    mediany = c(median(salariesFFSanF$TotalPay[salariesFFSanF$Year==2012]),
median(salariesFFSanF$TotalPay[salariesFFSanF$Year==2013]),
median(salariesFFSanF$TotalPay[salariesFFSanF$Year==2014])))

4.5 Wykres pokazujący zmiany średniej i mediany płac dla strażaków w San Francisco w latach 2012-2014

ggplot(sr_mediany_lata, aes(Rok, średnie)) + 
    geom_ribbon(aes(ymin = średnie, ymax = mediany),fill = "purple")+
  geom_line(aes(x=Rok,y=mediany,color="Mediany"),lwd=1)+
  geom_line(aes(x=Rok,y=średnie,color="Średnie"),lwd=1)+
  labs(x = "Year", y = "Total Pay (in Dollars)", title = "Zmiany różnicy między średnią i medianą na przestrzeni trzech lat wśród strażaków San Francisco")+
   theme_minimal()

Możemy zaobserwować zmiany w średniej i medianie wynagrodzeń strażaków w analizowanych latach, jak wspomniano wcześniej. Co więcej, różnica między medianą a średnią staje się coraz bardziej widoczna, przy czym średnia wypada na niższym poziomie w stosunku do mediany. To wskazuje, że większość strażaków otrzymuje wynagrodzenie wyższe niż przeciętna płaca w tej grupie zawodowej.

4.6 Jak płaca podstawowa, wynagrodzenie za nadgodziny i świadczenia są rozdzielane pomiędzy różne grupy?

salariesFFSanF$Benefits<-as.numeric(salariesFFSanF$Benefits)
salariesFFSanF$OvertimePay<-as.numeric(salariesFFSanF$OvertimePay)
salariesFFSanF$BasePay<-as.numeric(salariesFFSanF$BasePay)

4.7 Płaca podstawowa

sr_mediana_linie_basepay <- data.frame(nazwa = c("Średnia dla FT",
                                     "Mediana dla FT"), 
                                     wartość=c(mean(salariesFFSanF$BasePay[salariesFFSanF$Status=='FT']),
 median(salariesFFSanF$BasePay[salariesFFSanF$Status=='FT'])))

ggplot(salariesFFSanF,aes(x = BasePay)) +
  geom_histogram(binwidth=25000, color="green",fill="blue") +
  labs(x = "Base Pay", y = "Frequency", title = "Histogram płacy podstawowej dla strażaków w San Francisco") +
  geom_vline(data=sr_mediana_linie_basepay,aes(xintercept = wartość,
                               linetype = nazwa,
                               col = nazwa),linewidth=1,color="black")+
  theme_minimal()+
  facet_grid(. ~ Status)

## Warning: Removed 11 rows containing non-finite outside the scale range
## (`stat_bin()`).

4.8 Wynagrodzenie za nadgodziny

sr_mediana_linie_overtime <- data.frame(nazwa = c("Średnia",
                                     "Mediana"), wartość = c(mean(salariesFFSanF$OvertimePay),
 median(salariesFFSanF$OvertimePay)))

ggplot(salariesFFSanF,aes(x = OvertimePay)) +
  geom_histogram(binwidth=25000, color="green",fill="blue") +
  labs(x = "Overtime Pay", y = "Frequency", title = "Histogram płacy za nadgodziny dla strażaków w San Francisco") +
  geom_vline(data=sr_mediana_linie_overtime,aes(xintercept = wartość,
                               linetype = nazwa,
                               col = nazwa),linewidth=1,color="black")+
  theme_minimal()+
  facet_grid(. ~ Status)

4.9 Benefity

sr_mediana_linie_benefits <- data.frame(nazwa = c("Średnia",
                                     "Mediana"), wartość = c(mean(salariesFFSanF$Benefits),
 median(salariesFFSanF$Benefits)))

ggplot(salariesFFSanF,aes(x = Benefits)) +
  geom_histogram(binwidth=10000, color="red",fill="blue") +
  labs(x = "Benefits", y = "Frequency", title = "Histogram świadczeń/Benefitów dla strażaków w San Francisco") +
  geom_vline(data=sr_mediana_linie_benefits,aes(xintercept = wartość,
                               linetype = nazwa,
                               col = nazwa),linewidth=1,color="black")+
  theme_minimal()+
  facet_grid(. ~ Status)

Możemy zauważyć, że budżet przeznaczony na wynagrodzenie podstawowe, płatności za nadgodziny oraz świadczenia różni się w zależności od tego, czy strażacy są zatrudnieni w pełnym wymiarze godzin (FT), czy na pół etatu (ST). Zgodnie z przewidywaniami, wszystkie te wydatki są niższe w przypadku strażaków pracujących w niepełnym wymiarze godzin.

4.10 Czy w tym zestawie danych istnieją dowody na dyskryminację płacową ze względu na płeć?

Nie jesteśmy w stanie obliczyć tego dla tego zbioru danych, ponieważ nie mamy w nim zmiennej dotyczącej płci.

4.11 Oszacowanie gęstości rozkładu wynagrodzeń strażaków w San Francisco:

# Kernel density estimation
d <- density(salariesFFSanF$TotalPay,
             kernel = "gaussian",
             bw = 10000)
# Kernel density plot
plot(d, lwd = 2, main = "Gaussian kernel")

Dzięki zwiększeniu szerokości pasma do 10000,udało się wygładzić rozkład, który ułatwia nam wyciągnięcie wniosków. Mamy tutaj przykład rozkładu bimodalnego, który pomaga zrozumieć wcześniej zauważone zjawisko, że mediana przewyższa średnią. Wynika to z obecności znacznej grupy strażaków (stanowiącej drugie lokalne maksimum), którzy zarabiają niewielkie kwoty w przedziale <0; 50 000 USD>. Ich zarobki znacząco obniżają średnią płac, ale ich liczba nie jest na tyle duża, aby w istotny sposób wpłynąć na obniżenie mediany.