Kernel

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"

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)

Wykresy pudełkowe

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

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")

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:

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")

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")

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") 

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.

Ć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 #
#remotes::install_github("hericks/KDE") #tylko raz! potem #
#install.packages("htmltools")
library(KDE)
shiny_kde() 

Shiny applications not supported in static R Markdown documents

Tak, widać istotne różnice. Zwiększenie sd spłaszcza nam rozkład, modyfikacja średniej przesuwa środek rozkładu. Modyfikacja szerokości pasma wpływa na bliskość dopasowania rozkładu do danych (czym mniejsza wartość, tym bardziej dopasowana krzywa)

Ć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.

library(tidyverse)

## Warning: pakiet 'ggplot2' został zbudowany w wersji R 4.3.3

## Warning: pakiet 'readr' został zbudowany w wersji R 4.3.3

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

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

## [1] 2359   13

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="yellow",fill="red") +
  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="black")+
  theme_minimal()

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

##Poniżej histogram zmian wynagrodzeń strażaków z SF na przestrzeni lat, ponieważ działamy cały czas tylko na strażakach.

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="yellow",fill="red") +
  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="black")+
  theme_minimal()+
  facet_grid(. ~ Year)

Widzimy, że względem roku 2012, w 2013 roku liczba strażaków zarabiajacych tyle co srednia i mediana z 2012 roku zwiększyła się. Za to w roku 2014 wróciła do podobnych poziomów co w 2012 roku. Dowodem jest zmiana średniej płący na przestrzeni tych lat, gdzie po peaku w 2013 roku, spadła ona w 2014 roku do poziomu niższego niż w 2012

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])))

##Wykres pokazujący zmiany średniej i mediany płac dla strażaków w San Francisco na przestrzeni lat 2012-2014

ggplot(sr_mediany_lata, aes(Rok, średnie)) + 
    geom_ribbon(aes(ymin = średnie, ymax = mediany),fill = "yellow")+
  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")

Możemy zauwazyć zmiany wysokości średniej i mediany płacy dla strażaków na przestrzeni lat, tak jak było to wspomniane w poprzednim akapicie. Co więcej, widoczna jest zwiększająca się różnica między medianą a średnią, która jest coraz bardziej mniejsza od mediany. W związku z tym, możemy stwierdzić, że ponad połowa strażaków zarabia więcej niż średnia wynagrodzenia dla wszystkich strażaków.

• 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)

##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="orange",fill="darkred") +
  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()`).

##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="orange",fill="maroon") +
  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)

##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="darkorange") +
  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 wydatki budżetowe na płacę podstawową, płacę za nadgodziny i benefity różnią się w zależności od tego czy ktoś pracuje na pełen etat (FT) czy na pół etatu (ST). Zgodnie z oczekiwaniem, wszystkie te wartości są niższe dla strażaków, którzy nie pracują na pełen etat.

• Czy w tym zestawie danych istnieją dowody na dyskryminację płacową ze względu na płeć? W tym zestawie danych nie mamy informacji na temat płci, więc nie jesteśmy w stanie tego obliczyć.

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

# Kernel density estimation
d <- density(salariesFFSanF$TotalPay,
             kernel = "gaussian")

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

Wybrałem domyślne, gaussowskie jądro, które tak jak było wspomniane, zazwyczaj jest najlepsze. Mamy tu do czynienia z rozkładem bimodalnym, który niejako wyjaśnia nam zauważony wcześniej fakt, że mediana jest wyższa od średniej. Poprzez istotną liczbę (niejako drugą modę) strażaków zarabiających małe pieniądze <0;50 000 USD>, którzy zaniżają znacząco średnią wynagrodzeń, a nie jest ich wystarczająco dużo żeby istotnie zaniżyć medianę.