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:

# 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:

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ństwa moż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?

dane = TotalPay

# dane
dane<-salaries$TotalPay

# Kernel density estimation

gestosc <- density(dane)

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

# Kernel density estimation
gestosc_r <- density(dane,
             kernel = "rectangular")

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

# Kernel density estimation
gestosc_t <- density(dane,
             kernel = "triangular")

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

# Kernel density estimation
gestosc_e <- density(dane,
             kernel = "epanechnikov")

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

# Kernel density estimation
gestosc_b <- density(dane,
             kernel = "biweight")

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

dane = TotalPayBenefits

# Kernel density estimation
dane2<-salaries$TotalPayBenefits

gestosc2_e <- density(dane2,
             kernel = "epanechnikov")

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

# Kernel density estimation
gestosc2_r <- density(dane2,
             kernel = "rectangular")

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

# Kernel density estimation

gestosc2 <- density(dane2)

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

# Kernel density estimation
gestosc_nrd <- density(dane,
             bw = "nrd")

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

# Kernel density estimation
gestosc2_nrd <- density(dane2,
             bw = "nrd")

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

### Nieobciążona cross-walidacja

# Kernel density estimation
gestosc_ucv <- density(dane,
             bw = "ucv")

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

### Nieobciążona cross-walidacja

# Kernel density estimation
gestosc2_ucv <- density(dane2,
             bw = "ucv")
## Warning in bw.ucv(x): odnotowano wartość minimalną na jednym z końców zakresu
# Kernel density plot
plot(gestosc2_ucv, lwd = 2, main = "ucv bandwidth")

### Obciążona cross-walidacja

# Kernel density estimation
gestosc_bcv <- density(dane,
             bw = "bcv")

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

# Kernel density estimation
gestosc2_bcv <- density(dane2,
             bw = "bcv")

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

### Metoda Sheather & Jones (1991)
# Data

# Kernel density estimation
gestosc_sj <- density(dane,
             bw = "SJ")

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

### Metoda Sheather & Jones (1991)
# Data

# Kernel density estimation
gestosc2_sj <- density(dane2,
             bw = "SJ")

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

Wnioski:

Domyślna metoda daje najgładszą krzywą, która ignoruje drobne szczegóły, co może być zaletą w przypadku szumów w danych, ale wadą, jeśli chcemy zobaczyć szczegóły w rozkładzie. Metoda Silvermana wygładza krzywą bardziej niż Sheather-Jones, ale może być mniej trafna dla rozkładów nieregularnych. Metoda Sheather-Jones może wykazywać większą zmienność i wyraźniej ukazywać szczegóły, ale przy małych próbach może wprowadzać zbyt wiele fluktuacji, co może wyglądać jak szum.

Jądro prostokątne ma ostre krawędzie, co prowadzi do mniej płynnych oszacowań. Jądro Epanechnikova, biweight pokazują płynniejsze wykresy.

Funkcja gęstości staje się bardziej szczegółowa, co prowadzi do odwzorowania drobnych szczegółów w danych. W przypadku zbyt małej szerokości pasma, odwzorowuje szum zamiast rzeczywistego rozkładu. Występuje tendencja do wyłaniania wielu małych “górek”, które mogą nie mieć sensu w kontekście rzeczywistego rozkładu.

Duża szerokość pasma powoduje, że funkcja gęstości staje się bardziej wygładzona. Może to prowadzić do utraty drobnych szczegółów, takich jak lokalne maksima i minima. W przypadku bardzo dużej szerokości pasma można niedoszacować zmienności, co prowadzi do zbyt ogólnej estymacji, która nie odwzorowuje wszystkich cech rozkładu (np. wielomodalności).

Gdy dane są asymetryczne lub wielomodalne, metoda Sheather-Jones jest raczej najlepsza, ponieważ adaptuje się lepiej do nieregularności w danych. Zapewnia ona lepsze dopasowanie w przypadkach, gdy rozkład danych nie jest symetryczny. W porównaniu do bardziej wygładzonych metod, takich jak Silverman, Sheather-Jones lepiej radzi sobie z odwzorowaniem lokalnych max. i min.

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

# Filtracja danych dla 'JobTitle' i 'Agency'
firefighters <- subset(salaries, JobTitle == "FIREFIGHTER")
dane_strazacy<-firefighters$TotalPay
# Kernel density estimation

gestosc_f <- density(dane_strazacy)

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

 #podstawowy wykres KDE
ggplot(firefighters, aes(x = TotalPay)) +
  geom_density(fill = "blue", alpha = 0.3) +
  geom_histogram(aes(y = ..density..), bins = 20, fill = "grey", alpha = 0.5) +
  labs(title = "Kernel Density Estimation for Firefighters' Total Pay", x = "Firefighters' Total Pay", y = "Density")
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

###

# Przykład z różnymi jądrami oraz wykresem przekrojów

# gaussowski kernel
ggplot(firefighters, aes(x = TotalPay)) +
  geom_density(fill = "blue", alpha = 0.3, kernel = "gaussian", bw = 5000) +
  geom_histogram(aes(y = ..density..), bins = 20, fill = "grey", alpha = 0.5) +
  labs(title = "Kernel Density Estimation for Firefighters' Total Pay (Gaussian Kernel, bw=5000)", x = "Total Pay (in dollars)", y = "Density") +
  theme_minimal()

# epanechnikov
ggplot(firefighters, aes(x = TotalPay)) +
  geom_density(fill = "red", alpha = 0.3, kernel = "epanechnikov", bw = 10000) +
  geom_histogram(aes(y = ..density..), bins = 20, fill = "grey", alpha = 0.5) +
  labs(title = "Kernel Density Estimation for Firefighters' Total Pay (Epanechnikov Kernel, bw=10000)", x = "Total Pay (in dollars)", y = "Density") +
  theme_minimal()

ggplot(firefighters, aes(x = TotalPay)) +
  geom_density(fill = "blue", alpha = 0.3) +
  geom_histogram(aes(y = ..density..), bins = 20, fill = "grey", alpha = 0.5) +
  geom_vline(aes(xintercept = mean(TotalPay)), color = "red", linetype = "dashed", linewidth = 1) +
  geom_vline(aes(xintercept = median(TotalPay)), color = "green", linetype = "dashed", linewidth = 1) +
  labs(title = "Kernel Density Estimation with Mean and Median for Firefighters' Total Pay", x = "Total Pay (in dollars)", y = "Density") +
  theme_minimal()

Krzywa gęstości sugeruje, że wynagrodzenia strażaków mają skośny rozkład prawostronny. Istnieje kilka wyraźnych pików, wskazuje to na różne poziomy płac w zależności np. od rangi, doświadczenia lub dodatkowych godzin pracy.

Średnia zarobków znajduje się nieco na prawo od mediany. To pokazuje to, że istnieją wysokie wynagrodzenia, które podnoszą wartość średniej, podczas gdy większość danych (wynagrodzeń) koncentruje się w niższym zakresie.