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.
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ść.
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).
b <- density(runif(10))
f <- CDF(b)
f(0.5)
## [1] 0.6950212
plot(f)
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)
Przeczytaj artykuł naukowy “Kernel-smoothed cumulative distribution function estimation with akdensity” autorstwa Philippe Van Kerm.
Posłużymy się zbiorem danych diagnozy społecznej.
Na jego podstawie Twoim zadaniem jest oszacowanie rozkładu “p64 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")
library(ggplot2)
library(dplyr)
library(tidyr)
#install.packages("ggpubr")
library(ggpubr)
## Warning: pakiet 'ggpubr' został zbudowany w wersji R 4.4.2
##
## Dołączanie pakietu: 'ggpubr'
## Następujące obiekty zostały zakryte z 'package:spatstat.geom':
##
## border, rotate
# Zrobiłem to za pomocą pakietu ggpubr
Zbadanie struktury danych:
str(diagnoza)
## 'data.frame': 38461 obs. of 36 variables:
## $ imie_2011 : chr "WERONIKA" "ERNEST" "SYLWIA" "MARIOLA" ...
## $ waga_2013_osoby: num 0.278 0.278 0.278 0.278 0.475 ...
## $ lata_nauki_2013: num NA 11 16 16 7 12 12 12 12 9 ...
## $ wiek2013 : num 4 13 23 27 78 48 48 26 23 17 ...
## $ plec : Factor w/ 2 levels "mężczyzna","kobieta": 2 1 2 2 1 1 2 1 2 2 ...
## $ wojewodztwo : Factor w/ 17 levels "BD/ND/FALA","Dolnośląskie",..: 14 14 14 14 2 2 2 2 2 2 ...
## $ eduk4_2013 : Factor w/ 5 levels "BD/ND/FALA","podstawowe i niższe",..: NA 3 5 5 2 3 3 4 4 3 ...
## $ status9_2013 : Factor w/ 10 levels "BD/ND/FALA","pracownicy sekt. publicznego",..: 10 3 3 3 7 3 5 3 10 8 ...
## $ gp3 : Factor w/ 7 levels "WSPANIAŁE","UDANE",..: NA 2 3 5 3 2 2 NA 2 4 ...
## $ gp29 : Factor w/ 2 levels "PRZYJEMNOSCI",..: NA 2 1 2 2 2 2 NA 2 2 ...
## $ gp54_01 : Factor w/ 7 levels "ZDECYDOWANIE TAK ",..: NA 3 3 4 3 5 4 NA 4 5 ...
## $ gp54_02 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 1 1 1 2 1 1 NA 1 4 ...
## $ gp54_03 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 1 2 2 3 2 3 NA 1 3 ...
## $ gp54_04 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 4 5 4 5 5 5 NA 5 4 ...
## $ gp54_05 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 6 5 7 7 7 7 NA 6 7 ...
## $ gp54_06 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 5 3 5 7 2 3 NA 5 4 ...
## $ gp54_07 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 5 7 4 7 2 4 NA 7 7 ...
## $ gp54_08 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 5 5 6 5 6 3 NA 4 4 ...
## $ gp54_09 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 5 4 3 3 6 5 NA 4 3 ...
## $ gp54_10 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 2 1 3 2 2 1 NA 4 4 ...
## $ gp54_11 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 2 1 2 5 2 3 NA 4 4 ...
## $ gp54_12 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 1 5 3 4 6 5 NA 6 4 ...
## $ gp54_13 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 7 7 7 5 5 6 NA 6 3 ...
## $ gp54_14 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 1 1 1 3 2 2 NA 2 1 ...
## $ gp54_15 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 4 5 3 6 3 4 NA 6 4 ...
## $ gp54_16 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 1 1 1 2 4 2 NA 5 3 ...
## $ gp54_17 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 4 7 5 6 5 5 NA 7 3 ...
## $ gp54_18 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 1 1 1 3 3 4 NA 5 4 ...
## $ gp54_19 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 7 7 7 5 6 6 NA 6 4 ...
## $ gp54_20 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 1 3 2 1 6 3 NA 1 1 ...
## $ gp54_21 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 2 5 3 5 2 3 NA 4 1 ...
## $ gp54_22 : Factor w/ 7 levels "ZDECYDOWANIE TAK",..: NA 3 3 1 4 6 3 NA 1 5 ...
## $ gp60 : num NA 176 174 173 158 173 167 NA 167 175 ...
## $ gp61 : num NA 70 90 70 66 78 68 NA 64 65 ...
## $ gp64 : num NA 900 NA 950 2700 2500 1500 NA NA NA ...
## $ gp113 : num NA NA NA NA NA 40 15 NA NA NA ...
head(diagnoza)
## imie_2011 waga_2013_osoby lata_nauki_2013 wiek2013 plec wojewodztwo
## 1 WERONIKA 0.277653 NA 4 kobieta Świętokrzyskie
## 2 ERNEST 0.277653 11 13 mężczyzna Świętokrzyskie
## 3 SYLWIA 0.277653 16 23 kobieta Świętokrzyskie
## 4 MARIOLA 0.277653 16 27 kobieta Świętokrzyskie
## 5 WACŁAW 0.475255 7 78 mężczyzna Dolnośląskie
## 6 PIOTR 0.251773 12 48 mężczyzna Dolnośląskie
## eduk4_2013 status9_2013 gp3
## 1 <NA> inni bierni zawodowo <NA>
## 2 zasadnicze zawodowe/gimanzjum pracownicy sekt. prywatnego UDANE
## 3 wyższe i policealne pracownicy sekt. prywatnego DOSYĆ DOBRE
## 4 wyższe i policealne pracownicy sekt. prywatnego NIEZBYT UDANE
## 5 podstawowe i niższe emeryci DOSYĆ DOBRE
## 6 zasadnicze zawodowe/gimanzjum pracownicy sekt. prywatnego UDANE
## gp29 gp54_01 gp54_02 gp54_03
## 1 <NA> <NA> <NA> <NA>
## 2 POCZUCIE SENSU RACZEJ TAK ZDECYDOWANIE TAK ZDECYDOWANIE TAK
## 3 PRZYJEMNOSCI RACZEJ TAK ZDECYDOWANIE TAK TAK
## 4 POCZUCIE SENSU ANI TAK, ANI NIE ZDECYDOWANIE TAK TAK
## 5 POCZUCIE SENSU RACZEJ TAK TAK RACZEJ TAK
## 6 POCZUCIE SENSU RACZEJ NIE ZDECYDOWANIE TAK TAK
## gp54_04 gp54_05 gp54_06 gp54_07
## 1 <NA> <NA> <NA> <NA>
## 2 ANI TAK, ANI NIE NIE RACZEJ NIE RACZEJ NIE
## 3 RACZEJ NIE RACZEJ NIE RACZEJ TAK ZDECYDOWANIE NIE
## 4 ANI TAK, ANI NIE ZDECYDOWANIE NIE RACZEJ NIE ANI TAK, ANI NIE
## 5 RACZEJ NIE ZDECYDOWANIE NIE ZDECYDOWANIE NIE ZDECYDOWANIE NIE
## 6 RACZEJ NIE ZDECYDOWANIE NIE TAK TAK
## gp54_08 gp54_09 gp54_10 gp54_11
## 1 <NA> <NA> <NA> <NA>
## 2 RACZEJ NIE RACZEJ NIE TAK TAK
## 3 RACZEJ NIE ANI TAK, ANI NIE ZDECYDOWANIE TAK ZDECYDOWANIE TAK
## 4 NIE RACZEJ TAK RACZEJ TAK TAK
## 5 RACZEJ NIE RACZEJ TAK TAK RACZEJ NIE
## 6 NIE NIE TAK TAK
## gp54_12 gp54_13 gp54_14 gp54_15
## 1 <NA> <NA> <NA> <NA>
## 2 ZDECYDOWANIE TAK ZDECYDOWANIE NIE ZDECYDOWANIE TAK ANI TAK, ANI NIE
## 3 RACZEJ NIE ZDECYDOWANIE NIE ZDECYDOWANIE TAK RACZEJ NIE
## 4 RACZEJ TAK ZDECYDOWANIE NIE ZDECYDOWANIE TAK RACZEJ TAK
## 5 ANI TAK, ANI NIE RACZEJ NIE RACZEJ TAK NIE
## 6 NIE RACZEJ NIE TAK RACZEJ TAK
## gp54_16 gp54_17 gp54_18 gp54_19
## 1 <NA> <NA> <NA> <NA>
## 2 ZDECYDOWANIE TAK ANI TAK, ANI NIE ZDECYDOWANIE TAK ZDECYDOWANIE NIE
## 3 ZDECYDOWANIE TAK ZDECYDOWANIE NIE ZDECYDOWANIE TAK ZDECYDOWANIE NIE
## 4 ZDECYDOWANIE TAK RACZEJ NIE ZDECYDOWANIE TAK ZDECYDOWANIE NIE
## 5 TAK NIE RACZEJ TAK RACZEJ NIE
## 6 ANI TAK, ANI NIE RACZEJ NIE RACZEJ TAK NIE
## gp54_20 gp54_21 gp54_22 gp60 gp61 gp64 gp113
## 1 <NA> <NA> <NA> NA NA NA NA
## 2 ZDECYDOWANIE TAK TAK RACZEJ TAK 176 70 900 NA
## 3 RACZEJ TAK RACZEJ NIE RACZEJ TAK 174 90 NA NA
## 4 TAK RACZEJ TAK ZDECYDOWANIE TAK 173 70 950 NA
## 5 ZDECYDOWANIE TAK RACZEJ NIE ANI TAK, ANI NIE 158 66 2700 NA
## 6 NIE TAK NIE 173 78 2500 40
Wstępne przetwarzanie danych
# Usunięcie obserwacji z brakującymi wartościami w kluczowych kolumnach
diagnoza_clean <- diagnoza %>%
filter(!is.na(gp64), !is.na(wojewodztwo), !is.na(plec)) %>%
mutate(
plec = as.factor(plec), # Konwersja płci na typ faktorowy
wojewodztwo = as.factor(wojewodztwo) # Konwersja województwa na typ faktorowy
)
# Sprawdzenie rozmiaru i struktury przetworzonego zbioru danych
summary(diagnoza_clean)
## imie_2011 waga_2013_osoby lata_nauki_2013 wiek2013
## Length:19688 Min. : 0.0000 Min. : 0.00 Min. : 13.00
## Class :character 1st Qu.: 0.4816 1st Qu.:10.00 1st Qu.: 39.00
## Mode :character Median : 0.7919 Median :11.00 Median : 54.00
## Mean : 0.9970 Mean :11.79 Mean : 52.69
## 3rd Qu.: 1.2947 3rd Qu.:13.00 3rd Qu.: 65.00
## Max. :10.1024 Max. :28.00 Max. :101.00
## NA's :41 NA's :12
## plec wojewodztwo eduk4_2013
## mężczyzna: 8918 Mazowieckie :2489 BD/ND/FALA : 0
## kobieta :10770 Śląskie :1941 podstawowe i niższe :3949
## Wielkopolskie:1534 zasadnicze zawodowe/gimanzjum:5831
## Łódzkie :1449 średnie :5740
## Lubelskie :1415 wyższe i policealne :4116
## Dolnośląskie :1381 NA's : 52
## (Other) :9479
## status9_2013 gp3
## emeryci :6288 UDANE :7650
## pracownicy sekt. prywatnego :5177 DOSYĆ DOBRE :6949
## pracownicy sekt. publicznego:2666 ANI DOBRE, ANI ZŁE:3234
## renciści :1785 NIEZBYT UDANE : 992
## rolnicy :1263 WSPANIAŁE : 612
## (Other) :2433 (Other) : 224
## NA's : 76 NA's : 27
## gp29 gp54_01 gp54_02
## PRZYJEMNOSCI : 7389 NIE :5304 TAK :9084
## POCZUCIE SENSU:12221 ANI TAK, ANI NIE :4216 RACZEJ TAK :4143
## NA's : 78 RACZEJ NIE :2944 ZDECYDOWANIE TAK:3733
## RACZEJ TAK :2412 ANI TAK, ANI NIE:1772
## TAK :2374 NIE : 395
## (Other) :2368 (Other) : 484
## NA's : 70 NA's : 77
## gp54_03 gp54_04 gp54_05
## TAK :5093 NIE :4537 NIE :9390
## RACZEJ TAK :4609 ANI TAK, ANI NIE:4259 ZDECYDOWANIE NIE:3963
## ANI TAK, ANI NIE:3574 RACZEJ NIE :3541 RACZEJ NIE :3006
## ZDECYDOWANIE TAK:2305 RACZEJ TAK :3416 ANI TAK, ANI NIE:1460
## RACZEJ NIE :1918 TAK :2189 TAK : 795
## (Other) :2118 (Other) :1671 (Other) :1004
## NA's : 71 NA's : 75 NA's : 70
## gp54_06 gp54_07 gp54_08
## NIE :5161 NIE :5278 ANI TAK, ANI NIE:4720
## ANI TAK, ANI NIE:3464 ANI TAK, ANI NIE:4055 RACZEJ NIE :4030
## RACZEJ TAK :2917 ZDECYDOWANIE NIE:4015 RACZEJ TAK :3622
## RACZEJ NIE :2764 RACZEJ NIE :2806 NIE :3143
## TAK :2685 RACZEJ TAK :1542 TAK :2766
## (Other) :2628 (Other) :1904 (Other) :1330
## NA's : 69 NA's : 88 NA's : 77
## gp54_09 gp54_10 gp54_11
## ANI TAK, ANI NIE:4600 TAK :7293 TAK :5720
## TAK :4468 ZDECYDOWANIE TAK:4316 RACZEJ TAK :5600
## RACZEJ TAK :3689 RACZEJ TAK :3294 ANI TAK, ANI NIE:4171
## NIE :2571 ANI TAK, ANI NIE:2279 ZDECYDOWANIE TAK:1883
## ZDECYDOWANIE NIE:1521 NIE :1072 RACZEJ NIE :1072
## (Other) :2758 (Other) :1360 (Other) :1176
## NA's : 81 NA's : 74 NA's : 66
## gp54_12 gp54_13 gp54_14
## ANI TAK, ANI NIE:5530 TAK :4623 TAK :7480
## NIE :3789 ANI TAK, ANI NIE:3890 RACZEJ TAK :4383
## RACZEJ NIE :3326 NIE :3538 ZDECYDOWANIE TAK:4345
## TAK :2802 RACZEJ TAK :3159 ANI TAK, ANI NIE:2385
## RACZEJ TAK :2353 RACZEJ NIE :1975 NIE : 487
## (Other) :1804 (Other) :2429 (Other) : 540
## NA's : 84 NA's : 74 NA's : 68
## gp54_15 gp54_16 gp54_17
## ANI TAK, ANI NIE:4680 TAK :7490 ANI TAK, ANI NIE:4427
## NIE :4626 RACZEJ TAK :5158 TAK :4097
## RACZEJ NIE :4087 ANI TAK, ANI NIE:2641 NIE :3378
## RACZEJ TAK :3101 ZDECYDOWANIE TAK:2565 RACZEJ TAK :2934
## TAK :1852 RACZEJ NIE : 874 RACZEJ NIE :2522
## (Other) :1269 (Other) : 879 (Other) :2243
## NA's : 73 NA's : 81 NA's : 87
## gp54_18 gp54_19 gp54_20
## TAK :5380 NIE :5805 TAK :7708
## RACZEJ TAK :3832 ANI TAK, ANI NIE:3897 RACZEJ TAK :4265
## ANI TAK, ANI NIE:3131 RACZEJ NIE :2865 ZDECYDOWANIE TAK:3917
## ZDECYDOWANIE TAK:3091 ZDECYDOWANIE NIE:2701 ANI TAK, ANI NIE:2373
## RACZEJ NIE :1939 RACZEJ TAK :2062 RACZEJ NIE : 669
## (Other) :2234 (Other) :2270 (Other) : 683
## NA's : 81 NA's : 88 NA's : 73
## gp54_21 gp54_22 gp60
## ANI TAK, ANI NIE:4538 TAK :5633 Min. :110.0
## RACZEJ NIE :4347 ANI TAK, ANI NIE:4707 1st Qu.:162.0
## RACZEJ TAK :3851 RACZEJ TAK :4485 Median :168.0
## NIE :3388 ZDECYDOWANIE TAK:2177 Mean :168.6
## TAK :2068 RACZEJ NIE :1147 3rd Qu.:175.0
## (Other) :1420 (Other) :1399 Max. :210.0
## NA's : 76 NA's : 140 NA's :130
## gp61 gp64 gp113
## Min. : 30.00 Min. : 100 Min. : 0.00
## 1st Qu.: 65.00 1st Qu.: 1000 1st Qu.:40.00
## Median : 75.00 Median : 1500 Median :40.00
## Mean : 75.58 Mean : 1739 Mean :40.96
## 3rd Qu.: 85.00 3rd Qu.: 2000 3rd Qu.:45.00
## Max. :180.00 Max. :50000 Max. :90.00
## NA's :84 NA's :9272
Estymacja rozkładu gęstości dla dochodu według województw i płci:
# Gęstość dochodu wg płci i województw
ggplot(diagnoza_clean, aes(x = gp64, color = plec, fill = plec)) +
geom_density(alpha = 0.4) +
facet_wrap(~ wojewodztwo) +
labs(
title = "Rozkład gęstości dochodu miesięcznego netto wg województwa i płci",
x = "Dochód netto (gp64)",
y = "Gęstość"
) +
theme_minimal()
# Dystrybuanta dochodu wg płci i województw
ggplot(diagnoza_clean, aes(x = gp64, color = plec, group = plec)) +
stat_ecdf(geom = "step") +
facet_wrap(~ wojewodztwo) +
labs(
title = "Dystrybuanta dochodu miesięcznego netto wg województwa i płci",
x = "Dochód netto (gp64)",
y = "Prawdopodobieństwo skumulowane"
) +
theme_minimal()
Podsumowanie danych:
# Statystyki opisowe
summary_stats <- diagnoza_clean %>%
group_by(wojewodztwo, plec) %>%
summarise(
srednia_dochod = mean(gp64, na.rm = TRUE),
mediana_dochod = median(gp64, na.rm = TRUE),
odchylenie_std = sd(gp64, na.rm = TRUE),
liczba_obserwacji = n(),
.groups = "drop"
)
print(summary_stats)
## # A tibble: 32 × 6
## wojewodztwo plec srednia_dochod mediana_dochod odchylenie_std
## <fct> <fct> <dbl> <dbl> <dbl>
## 1 Dolnośląskie mężczyzna 2130. 1800 1358.
## 2 Dolnośląskie kobieta 1627. 1400 992.
## 3 Kujawsko-Pomorskie mężczyzna 1887. 1500 2469.
## 4 Kujawsko-Pomorskie kobieta 1327. 1200 725.
## 5 Lubelskie mężczyzna 1761. 1500 1396.
## 6 Lubelskie kobieta 1325. 1100 774.
## 7 Lubuskie mężczyzna 1986. 1800 1251.
## 8 Lubuskie kobieta 1541. 1400 832.
## 9 Łódzkie mężczyzna 1936. 1500 1736.
## 10 Łódzkie kobieta 1447. 1200 988.
## # ℹ 22 more rows
## # ℹ 1 more variable: liczba_obserwacji <int>
Najwyższy średni dochód w przypadku mężczyzn wystąpił w województwie pomorskim (2279), a w przypadku kobiet w województwie mazowieckim (1780). Natomiast najniższy dochód wystąpił w zarówno w przypadku mężczyzn (1588) jak i w przypadku kobiet (1239) w woj. podkarpackim.