Kernel

## 
## Dołączanie pakietu: 'dplyr'

## Następujące obiekty zostały zakryte z 'package:stats':
## 
##     filter, lag

## Następujące obiekty zostały zakryte z 'package:base':
## 
##     intersect, setdiff, setequal, union

## Ładowanie wymaganego pakietu: viridisLite

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?

3.1 Wnioski

Zauważone podczas eksperymentowania z różnymi parametrami estymatora funkcji gęstości kernel:

Rozkłady:

Rozkład normalny ma dzwonowaty kształt, skoncentrowany wokół średniej. Zmiana średniej przesuwa wykres w prawo lub lewo, a większe odchylenie standardowe spłaszcza go.

Rozkład beta może przybierać różne, często asymetryczne kształty, zależne od parametrów α i β.

Rozkład wykładniczy może sprawiać że gęstości szybko maleje od początkowej wartości. Im większy parametr, tym szybciej spada.

Rozkład jednostajny przypisuje jednakowe prawdopodobieństwo dla wszystkich wartości w określonym przedziale.

Nietypowe rozkłady mają nieregularny kształt i często zależą od specyficznych danych.

Kernele

Rodzaj decyduje o sposobie wygładzania: Gaussowskie daje dzwonowaty,gładki kształt. Biweight więcej uwagi przypisuje punktom blisko środka. Cosinusowe wygładza wartości stopniowo w kierunku krańców przedziału.

Szerokość pasma wpływa na precyzję i gładkość wykresu:

Węższe pasmo umożliwia lepsze odwzorowanie detali, ale może prowadzić do zbyt szczegółowego dopasowania. Optymalne pasmo dobrze balansuje szczegóły i gładkość. Szersze pasmo daje gładszy wykres, ale może zgubić niektóre detale danych.

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

Shiny applications not supported in static R Markdown documents

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.

4.1 Poprzednie biblioteczki

4.2 Badanie struktury danych Wynagrodzeń (zbiór danych salaries) strażaków w San Francisco

salaries3a<- subset(salaries, JobTitle=="Firefighter")
dim(salaries3a)

## [1] 2359   13

hist(salaries3a$TotalPay,main="Całkowite wynagrodzenie", xlab="Wynagrodzenia w dolarach",ylab="Częstotliwość")
abline(v = mean(salaries3a$TotalPay),lty="dashed")
abline(v = median(salaries3a$TotalPay))
legend("topright", legend=c("Mediana","Średnia"),lty=c("solid","dashed"))

salaries3a <- subset(salaries, JobTitle=="Firefighter")
dim(salaries3a)

## [1] 2359   13

Szukamy czy nie ma nigdzie NA

##   na_count
## 1        0

4.2.1 Szukamy czy nie ma wartości odstajacych

Wykres pudełkowy wynagrodzenia strażaków w San Francisco z wartościami odstającymi

ggplot(salaries3a, aes(y = TotalPay)) +
  geom_boxplot() +
  coord_cartesian(ylim = c(0, 300000)) 

Wykres pudełkowy wynagrodzenia strażaków w San Francisco besz wartości odstających

ggplot(salaries3a, aes(y = TotalPay)) +
  geom_boxplot(outlier.shape = NA) +
  coord_cartesian(ylim = c(0, 250000)) 

salaries3a %>% 
  filter(TotalPay < 250000 , TotalPay >75000 ) %>% 
  ggplot(aes(y = TotalPay)) +
  geom_boxplot() +
  coord_cartesian(ylim = c(0, 250000))

Tworzymy zbiór danych bez wartości odstających

salaries3<-salaries3a %>%
  filter(TotalPay < 250000 , TotalPay >75000 )

Aby sprawdzaić czy zbiór salaries 3 pasuje ona do naszej analizy , przejrzymy czy na jego podstawie da się odpowiedzieć an pytanie: jak zmieniały się wynagrodzenia wsród strazaków?

Wykres ukazujący zmian wynagrodzenia wsród strazaków na podstawie zbiory salaries3

ggplot(data=salaries3, aes(x=BasePay, group=Year, fill=Year)) +
    geom_density(adjust=1.5) +
    theme_ipsum() +
    facet_wrap(~Year) +
    theme(
      legend.position="none",
      panel.spacing = unit(0.1, "lines"),
      axis.ticks.x=element_blank()
    )

Wykres ukazujący zmian wynagrodzenia wsród strazaków na podstawie zbiory salaries3a

ggplot(data=salaries3a, aes(x=BasePay, group=Year, fill=Year)) +
    geom_density(adjust=1.5) +
    theme_ipsum() +
    facet_wrap(~Year) +
    theme(
      legend.position="none",
      panel.spacing = unit(0.1, "lines"),
      axis.ticks.x=element_blank()
    )

Decydujemy się pozostać przy zbiore salaries3a aby ukazać bardziej zróżnicowe wyniki w poszczgólnych latach.

4.3 Estymacji funkcji gęstości dla wynagrodzeń strażaków w San Francisco

Histogram i estymacja funkcji gęstości dla Wynagrodzenie całkowitego strażaków

ggplot(salaries3a, aes(x = TotalPay)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30
                 , fill = "green", color = "black", alpha = 0.2) +
  geom_density(color = "purple", size = 1) +
  labs(#title = "Histogram i estymacja funkcji gęstości dla Wynagrodzenie całkowitego strażaków",
       x = "Wynagrodzenie całkowite strażaków",
       y = "Gęstość") +
  theme_minimal()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

4.4 Estymacji funkcji gęstości dla wynagrodzeń strażaków w San Francisco za pomoca kernela

data1 <- salaries3a$TotalPay
d1 <- density(data1,
             kernel = "cosine")
plot(d1, lwd = 2, main = "Cosine kernel")

d2 <- density(data1,
             bw = "SJ")
plot(d2, lwd = 2, main = "SJ bandwidth")

d3 <- density(data1,
             bw = "nrd")
plot(d3, lwd = 2, main = "nrd bandwidth")

4.4.1 Jak zmieniały się wynagrodzenia w czasie wsród strażaków?

Wykres wynagrodzenia za w zależności od roku

ggplot(data=salaries3a, aes(x=TotalPay, group=Year, fill=Year)) +
    geom_density(adjust=1.5, alpha=.4) +
    theme_ipsum()

4.4.2 Jak zmieniały się wynagrodzenia za nadgodziny w czasie wsród strażaków?

Wykres wynagrodzenia za nadgodziny w zależności od roku

ggplot(data=salaries3a, aes(x=OvertimePay, group=as.factor(Year), fill=Year)) +
    geom_density(adjust=1.5, alpha=0.4) +
    theme_ipsum()

Wykres wynagrodzenia podstawowego w zależności od roku

ggplot(data=salaries3, aes(x=BasePay, group=Year, fill=Year)) +
    geom_density(adjust=1.5) +
    theme_ipsum() +
    facet_wrap(~Year) +
    theme(
      legend.position="none",
      panel.spacing = unit(0.1, "lines"),
      axis.ticks.x=element_blank()
    )

grupy <- salaries3a %>%
  filter(Year %in% c("2012", "2013", "2014"))

# Tworzenie wykresu z poprawionym kolorowaniem
ggplot(grupy, aes(x = TotalPay, color = as.factor(Year))) +
  geom_density(size = 1) +
  labs(title = "Estymacja funkcji gęstości dla wynagrodzeń w różnych latach wsród strażaków",
       x = "Wynagrodzenie (Total Pay)",
       y = "Gęstość",
       color = "ROK") +
  theme_minimal()

Wykres ukazujący zmian w dodatkach do wynagrodzenia wsród strazaków na podstawie zbiory salaries3

ggplot(data=salaries3a, aes(x=Benefits, group=Year, fill=Year)) +
    geom_density(adjust=1.5) +
    theme_ipsum() +
    facet_wrap(~Year) +
    theme(
      legend.position="none",
      panel.spacing = unit(0.1, "lines"),
      axis.ticks.x=element_blank()
    )

5 Pytania badawcze

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

salaries$JobTitle %>% table() %>% sort(decreasing = TRUE) %>%
head(20)

## .
##             Transit Operator                Special Nurse 
##                         7036                         4389 
##             Registered Nurse Public Svc Aide-Public Works 
##                         3736                         2518 
##             Police Officer 3                    Custodian 
##                         2421                         2418 
##             TRANSIT OPERATOR                  Firefighter 
##                         2388                         2359 
##            Recreation Leader       Patient Care Assistant 
##                         1971                         1945 
##               Deputy Sheriff               Police Officer 
##                         1933                         1476 
##                SPECIAL NURSE       Public Service Trainee 
##                         1402                         1328 
##             REGISTERED NURSE             Police Officer 2 
##                         1219                         1141 
##    Attorney (Civil/Criminal)                       Porter 
##                         1126                         1095 
##                   Sergeant 3              General Laborer 
##                         1047                         1033

grupy1 <- salaries %>%
  filter(JobTitle %in% c("Firefighter", "Police Officer", "Engineer", "Special Nurse","Custodian","NURSE","General Laborer"))


ggplot(grupy1, aes(x = TotalPay, color = JobTitle)) +
  geom_density(size = 1) +
  labs(title = "Estymacja funkcji gęstości dla wynagrodzeń  wsród różnych zawodów",
       x = "Wynagrodzenie",
       y = "Gęstość",
       color = "Zawód") +
  theme_minimal()