Nieklasyczne metody statystyki

Nieparametryczna estymacja dystrybuanty

1 Wprowadzenie

Cel: wyznaczenie obszaru ufności dla dystrybuanty nieznanego rozkładu, a nie tylko oszacowania parametrów, od jakich zależą jej wartości.

• Brzegi tego obszaru są wykresami funkcji „przedziałami stałych” (funkcji schodkowych).

• Jeżeli przy wyznaczaniu pasma ufności dla dystrybuanty otrzymamy lewy kraniec przedziału będący liczbą ujemną, to zastępujemy ją przez zero.

• Jeżeli otrzymamy prawy kraniec przedziału większy od jedności, to przyjmujemy, że jest on równy jeden.

• Określenie obszaru ufności dla dystrybuanty w przedstawiony sposób polega na wyznaczeniu przedziałowego oszacowania dla każdej wartości dystrybuanty.

1.1 Funkcja CDF

Funkcja w programie R odpowiedzialna za estymację to np. CDF z pakietu spatstat. CDF jest metodą ogólną, z metodą dla klasy “gęstość”.

Oblicza ona skumulowaną funkcję rozkładu, której gęstość prawdopodobieństwa została oszacowana i zapisana w obiekcie f. Obiekt f musi należeć do klasy “gęstość” i zazwyczaj zostałby uzyskany z wywołania funkcji gęstość.

1.2 Funkcja kde

Pakiet R o nazwie snpar zawiera kilka uzupełniających metod statystyki nieparametrycznej, w tym test kwantylowy, test trendu Coxa-Stuarta, test przebiegów, test normalnego wyniku, estymację jądra PDF i CDF, estymację regresji jądra i test jądra Kołmogorowa-Smirnowa.

Funkcja kde zawiera obliczanie zarówno nieparametrycznego estymatora jądra funkcji gęstości prawdopodobieństwa (PDF) jak i funkcji rozkładu skumulowanego (CDF).

1.3 Przykład 1.

   b <- density(runif(10))
   f <- CDF(b)
   f(0.5)

## [1] 0.4280671

   f(1)

## [1] 0.9278786

   plot(f)

1.4 Przykład 2.

x <- rnorm(200,2,3)
# with default bandwidth
# kde(x, kernel = "quar", plot = TRUE)

# with specified bandwidth
# kde(x, h = 4, kernel = "quar", plot = TRUE)

1.5 Przeczytaj

Przeczytaj artykuł naukowy “Kernel-smoothed cumulative distribution function estimation with akdensity” autorstwa Philippe Van Kerm.

1.6 Zadanie

Posłużymy się zbiorem danych diagnozy społecznej.

Na jego podstawie Twoim zadaniem jest oszacowanie rozkładu “gp64 Pana/Pani wlasny (osobisty) dochod miesieczny netto (na reke)” według województw/płci.

Postaraj się oszacować zarówno rozkład gęstości jak i skumulowanej gęstości (dystrybuanty).

data("diagnoza")
data("diagnozaDict")

# Filtrujemy dane 
data_filtered <- diagnoza %>%
  select(gp64, wojewodztwo, plec) %>%
  filter(!is.na(gp64) & gp64 > 0 & wojewodztwo != "BD/ND/FALA") %>%
  mutate(wojewodztwo = as.factor(wojewodztwo),  # Województwo
         plec = as.factor(plec))  # Płeć

# Podział województw na dwie grupy
woj_group_1 <- levels(data_filtered$wojewodztwo)[1:4]
woj_group_2 <- levels(data_filtered$wojewodztwo)[5:8]
woj_group_3 <- levels(data_filtered$wojewodztwo)[9:12]
woj_group_4 <- levels(data_filtered$wojewodztwo)[13:16]

# Obliczanie rozkładów KDE i dystrybuant
# Estymacja dla każdej kombinacji województwa i płci
kde_results <- data_filtered %>%
  group_by(wojewodztwo, plec) %>%
  summarise(density = list(density(gp64, na.rm = TRUE)), .groups = "drop")

cdf_results <- data_filtered %>%
  group_by(wojewodztwo, plec) %>%
  summarise(
    cdf = list(ecdf(gp64)),
    .groups = "drop"
  )

1.7 Wykresy gęstości i dystrubuanty skumulowa z podziałem na 2 grupy województw (dla większej czytelności)

# Wykres rozkładu gęstości dla grupy 1
plot_density_group1 <- ggplot(data_filtered %>% filter(wojewodztwo %in% woj_group_1), aes(x = gp64, color = plec)) +
  geom_density() +
  facet_wrap(~ wojewodztwo, scales = "free_y") +
  labs(
    title = "Rozkład gęstości dochodu netto (gp64) - Grupa 1 województw",
    x = "Dochód netto (zł)",
    y = "Gęstość",
    color = "Płeć"
  ) +
  scale_x_continuous(labels = scales::comma, limits = c(0, quantile(data_filtered$gp64, 0.99))) +
  theme_minimal() +
  theme(plot.title = element_text(size = 16),
        axis.text = element_text(size = 12),
        axis.title = element_text(size = 14),
        strip.text = element_text(size = 14))

# Wykres rozkładu gęstości dla grupy 2
plot_density_group2 <- ggplot(data_filtered %>% filter(wojewodztwo %in% woj_group_2), aes(x = gp64, color = plec)) +
  geom_density() +
  facet_wrap(~ wojewodztwo, scales = "free_y") +
  labs(
    title = "Rozkład gęstości dochodu netto (gp64) - Grupa 2 województw",
    x = "Dochód netto (zł)",
    y = "Gęstość",
    color = "Płeć"
  ) +
  scale_x_continuous(labels = scales::comma, limits = c(0, quantile(data_filtered$gp64, 0.99))) +
  theme_minimal() +
  theme(plot.title = element_text(size = 16),
        axis.text = element_text(size = 12),
        axis.title = element_text(size = 14),
        strip.text = element_text(size = 14))
# Wykres rozkładu gęstości dla grupy 3
plot_density_group3 <- ggplot(data_filtered %>% filter(wojewodztwo %in% woj_group_3), aes(x = gp64, color = plec)) +
  geom_density() +
  facet_wrap(~ wojewodztwo, scales = "free_y") +
  labs(
    title = "Rozkład gęstości dochodu netto (gp64) - Grupa 3 województw",
    x = "Dochód netto (zł)",
    y = "Gęstość",
    color = "Płeć"
  ) +
  scale_x_continuous(labels = scales::comma, limits = c(0, quantile(data_filtered$gp64, 0.99))) +
  theme_minimal() +
  theme(plot.title = element_text(size = 16),
        axis.text = element_text(size = 12),
        axis.title = element_text(size = 14),
        strip.text = element_text(size = 14))
# Wykres rozkładu gęstości dla grupy 4
plot_density_group4 <- ggplot(data_filtered %>% filter(wojewodztwo %in% woj_group_4), aes(x = gp64, color = plec)) +
  geom_density() +
  facet_wrap(~ wojewodztwo, scales = "free_y") +
  labs(
    title = "Rozkład gęstości dochodu netto (gp64) - Grupa 4 województw",
    x = "Dochód netto (zł)",
    y = "Gęstość",
    color = "Płeć"
  ) +
  scale_x_continuous(labels = scales::comma, limits = c(0, quantile(data_filtered$gp64, 0.99))) +
  theme_minimal() +
  theme(plot.title = element_text(size = 16),
        axis.text = element_text(size = 12),
        axis.title = element_text(size = 14),
        strip.text = element_text(size = 14))


# Wykres dystrybuanty skumulowanej dla grupy 1
plot_cdf_group1 <- ggplot(data_filtered %>% filter(wojewodztwo %in% woj_group_1), aes(x = gp64, color = plec)) +
  stat_ecdf(geom = "step") +
  facet_wrap(~ wojewodztwo, scales = "free_y") +
  labs(
    title = "Dystrybuanta skumulowana dochodu netto (gp64) - Grupa 1 województw",
    x = "Dochód netto (zł)",
    y = "Prawdopodobieństwo skumulowane",
    color = "Płeć"
  ) +
  scale_x_continuous(labels = scales::comma, limits = c(0, quantile(data_filtered$gp64, 0.99))) +
  theme_minimal() +
  theme(plot.title = element_text(size = 16),
        axis.text = element_text(size = 12),
        axis.title = element_text(size = 14),
        strip.text = element_text(size = 14))

# Wykres dystrybuanty skumulowanej dla grupy 2
plot_cdf_group2 <- ggplot(data_filtered %>% filter(wojewodztwo %in% woj_group_2), aes(x = gp64, color = plec)) +
  stat_ecdf(geom = "step") +
  facet_wrap(~ wojewodztwo, scales = "free_y") +
  labs(
    title = "Dystrybuanta skumulowana dochodu netto (gp64) - Grupa 2 województw",
    x = "Dochód netto (zł)",
    y = "Prawdopodobieństwo skumulowane",
    color = "Płeć"
  ) +
  scale_x_continuous(labels = scales::comma, limits = c(0, quantile(data_filtered$gp64, 0.99))) +
  theme_minimal() +
  theme(plot.title = element_text(size = 16),
        axis.text = element_text(size = 12),
        axis.title = element_text(size = 14),
        strip.text = element_text(size = 14))
# Wykres dystrybuanty skumulowanej dla grupy 3
plot_cdf_group3 <- ggplot(data_filtered %>% filter(wojewodztwo %in% woj_group_3), aes(x = gp64, color = plec)) +
  stat_ecdf(geom = "step") +
  facet_wrap(~ wojewodztwo, scales = "free_y") +
  labs(
    title = "Dystrybuanta skumulowana dochodu netto (gp64) - Grupa 3 województw",
    x = "Dochód netto (zł)",
    y = "Prawdopodobieństwo skumulowane",
    color = "Płeć"
  ) +
  scale_x_continuous(labels = scales::comma, limits = c(0, quantile(data_filtered$gp64, 0.99))) +
  theme_minimal() +
  theme(plot.title = element_text(size = 16),
        axis.text = element_text(size = 12),
        axis.title = element_text(size = 14),
        strip.text = element_text(size = 14))
# Wykres dystrybuanty skumulowanej dla grupy 4
plot_cdf_group4 <- ggplot(data_filtered %>% filter(wojewodztwo %in% woj_group_4), aes(x = gp64, color = plec)) +
  stat_ecdf(geom = "step") +
  facet_wrap(~ wojewodztwo, scales = "free_y") +
  labs(
    title = "Dystrybuanta skumulowana dochodu netto (gp64) - Grupa 4 województw",
    x = "Dochód netto (zł)",
    y = "Prawdopodobieństwo skumulowane",
    color = "Płeć"
  ) +
  scale_x_continuous(labels = scales::comma, limits = c(0, quantile(data_filtered$gp64, 0.99))) +
  theme_minimal() +
  theme(plot.title = element_text(size = 16),
        axis.text = element_text(size = 12),
        axis.title = element_text(size = 14),
        strip.text = element_text(size = 14))

# Wyświetlanie wykresów
print(plot_density_group1)

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

print(plot_density_group2)

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

print(plot_density_group3)

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

print(plot_density_group4)

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

print(plot_cdf_group1)

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

print(plot_cdf_group2)

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

print(plot_cdf_group3)

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

print(plot_cdf_group4)

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

Analizując zarówno gęstości, jak i dystrybuanty rozkładów, zauważamy, że kobiety mają niższy dochód netto niż mężczyźni. Szczyty (mody) rozkładów gęstości dla kobiet znajdują się na niższych wartościach niż u mężczyzn. Podobne obserwacje możemy poczynić na wykresach dystrybuanty skumulowanej, gdzie krzywe dla kobiet rosną szybciej i poziomują się, osiągając 75% w okolicach 3 000 zł, podczas gdy u mężczyzn wartość dystrybuanty na tym samym poziomie wynosi około 50%. Różnice te zmieniają się w zależności od województwa, ale ogólnie można stwierdzić, że w każdym z nich mężczyźni mają wyższy dochód netto miesięczny.