Nieklasyczne metody statystyki

Nieparametryczna estymacja dystrybuanty

Zadanie

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

new_df <- data.frame(p64_var=diagnoza$gp64,
                     region=diagnoza$wojewodztwo,
                     gender=diagnoza$plec)

num_nan_values <- colSums(is.na(new_df))
for(col_name in names(num_nan_values)){
  cat(paste('Liczba braków danych w kolumnie', col_name, 'wynosi', num_nan_values[col_name], '\n'))
}

## Liczba braków danych w kolumnie p64_var wynosi 18773 
## Liczba braków danych w kolumnie region wynosi 0 
## Liczba braków danych w kolumnie gender wynosi 40

Liczba braków danych w kolumnie p64_var stanowi blisko połowę wszystkich danych tej zmiennej. Patrząc jedynie na nowy zbiór danych wyselekcjonowany z “Diagnozy” nie ma podstaw do określenia natury braków danych, ponieważ w przypadku pozostałych dwóch zmiennych braki danych stanowią nieznaczną część lub w ogóle ich nie ma (zmienna województwo). Nie mając pełnej informacji o rozkładzie zmiennej “p64 Pana/Pani wlasny (osobisty) dochod miesieczny netto (na reke)” należy pamiętać o uniknięciu generalizacji wniosków. Z tego powodu, wykonam dwie różne estymację. Pierwszą po usunięciu braków danych oraz drugą dopiero po przekształceniu braków danych na średnią wartość ze wszystkich dostępnych obserwacji.

new_df_clean <- na.omit(new_df)

cat(paste('Unikalne wartości dotyczące zmiennej płeć to:', unique(new_df_clean$gender), '\n'))

## Unikalne wartości dotyczące zmiennej płeć to: mężczyzna 
##  Unikalne wartości dotyczące zmiennej płeć to: kobieta

female_data <- new_df_clean$p64_var[new_df_clean$gender == "kobieta"]
male_data <- new_df_clean$p64_var[new_df_clean$gender == "mężczyzna"]

female_density <- density(female_data, bw=3000)
male_density <- density(male_data, bw=3000)

plot(female_density, lwd=2, main=paste('Wykres gęstości dla dochodów względem płci - kobieta')) 

plot(male_density, lwd=2, main=paste('Wykres gęstości dla dochodów względem płci - mężczyzna')) 

ggplot(new_df_clean, aes(x=p64_var, color=gender, fill=gender)) +
  geom_density(alpha=0.3, bw=4000) + 
  labs(title='Porównanie gęstości dochodów względem płci (usunięte nan)', 
       x='Zarobki',
       y='Gęstość')+
  theme_minimal()

female_cdf <- approxfun(female_density$x, cumsum(female_density$y) / sum(female_density$y))
male_cdf <- approxfun(male_density$x, cumsum(male_density$y) / sum(male_density$y))

cat(paste('Jaki procent respondentów będących kobietą osiąga dochód netto niższy lub równy 10000: ', female_cdf(10000)*100, '%', sep="", '\n'))

## Jaki procent respondentów będących kobietą osiąga dochód netto niższy lub równy 10000: 99.4422417105648%

cat(paste('Jaki procent respondentów będących mężczyzną osiąga dochód netto niższy lub równy 10000: ',male_cdf(10000)*100, '%', sep=""))

## Jaki procent respondentów będących mężczyzną osiąga dochód netto niższy lub równy 10000: 98.8422799360035%

plot(female_density$x, female_cdf(female_density$x), type = "l", col = "red",
     xlab = "Value", ylab = "CDF", main = "CDF for Female data")

plot(male_density$x, male_cdf(male_density$x), type = "l", col = "blue",
     xlab = "Value", ylab = "CDF", main = "CDF for Male data")

cat(paste('Częstości występowania wartości w zmiennej płeć:\n'))

## Częstości występowania wartości w zmiennej płeć:

table(new_df$gender)

## 
## mężczyzna   kobieta 
##     18355     20066

cat('\n')            

new_df_impute <- new_df %>%
  mutate(p64_var=ifelse(is.na(p64_var), mean(p64_var, na.rm=TRUE), p64_var),
         gender=ifelse(gender=='kobieta', 1, ifelse(gender=='mężczyzna', 2, NA)) %>%
           as.numeric(),
         gender=ifelse(is.na(gender), 1, gender))

cat(paste('Unikalne wartości dotyczące zmiennej płeć to:', unique(new_df_impute$gender), '\n'))

## Unikalne wartości dotyczące zmiennej płeć to: 1 
##  Unikalne wartości dotyczące zmiennej płeć to: 2

female_data_impute <- new_df_impute$p64_var[new_df_impute$gender == 1]
male_data_impute <- new_df_impute$p64_var[new_df_impute$gender == 2]

female_density_impute <- density(female_data_impute, bw=3000)
male_density_impute <- density(male_data_impute, bw=3000)

plot(female_density_impute, lwd=2, main=paste('Wykres gęstości dla dochodów względem płci - kobieta')) 

plot(male_density_impute, lwd=2, main=paste('Wykres gęstości dla dochodów względem płci - mężczyzna')) 

ggplot(new_df_impute, aes(x=p64_var, color=factor(gender, labels=c('mężczyzna', 'kobieta')), fill=factor(gender, labels=c('mężczyzna', 'kobieta')))) +
  geom_density(alpha=0.3, bw=3000) + 
  labs(title='Porównanie gęstości zarobków względem płci (imputacja)', 
       x='Zarobki',
       y='Gęstość',
       color='Płeć',
       fill='Płeć')+
  theme_minimal()

female_cdf_impute <- approxfun(female_density_impute$x, cumsum(female_density_impute$y) / sum(female_density_impute$y))
male_cdf_impute <- approxfun(male_density_impute$x, cumsum(male_density_impute$y) / sum(male_density_impute$y))

paste('Jaki procent respondentów będących kobietą osiąga dochód netto niższy lub równy 10000:', female_cdf_impute(10000)*100, '%')

## [1] "Jaki procent respondentów będących kobietą osiąga dochód netto niższy lub równy 10000: 99.5692651140786 %"

paste('Jaki procent respondentów będących mężczyzną osiąga dochód netto niższy lub równy 10000:',male_cdf_impute(10000)*100, '%')

## [1] "Jaki procent respondentów będących mężczyzną osiąga dochód netto niższy lub równy 10000: 99.2940367645413 %"

plot(female_density_impute$x, female_cdf_impute(female_density_impute$x), type = "l", col = "red",
     xlab = "Value", ylab = "CDF", main = "CDF for Female data")

plot(male_density_impute$x, male_cdf_impute(male_density_impute$x), type = "l", col = "blue",
     xlab = "Value", ylab = "CDF", main = "CDF for Male data")

Powyżej zaprezentowałem oszacowania rozkładów gęstości oraz skumulowanej gęstości z podziałem na płeć. Wykres, który powstał w wyniku usuwania brakujących wartości lepiej obrazuje różnice w dochodzie miesięcznym netto. Mniej widoczne jest to w przypadku danych, gdzie użyta została technika uzupełniania braków, co jest jak najbardziej prawdziwe, gdyż połowa braków została uzupełniona średnią (stąd m.in. mediana taka sama dla dwóch grup).

Korzystając z wykresu dystrybuanty można napisać, że : - 99.6% procent respondentów (kobiet) osiąga dochód netto niższy lub równy 10000 (dane po usunięciu NAN), - 99.4% procent respondentów (kobiet) osiąga dochód netto niższy lub równy 10000 (dane imputowane), - 99,3% procent respondentów (mężczyzn) osiąga dochód netto niższy lub równy 10000 (dane po usunięciu NAN), - 98,8% procent respondentów (mężczyzn) osiąga dochód netto niższy lub równy 10000 (dane imputowane).

Jak to wygląda dla województw? Wezmę województwa, z których pochodziło najwięcej ankietowanych tj. Mazowieckie, Śląśkie i Wielkopolskie. Przygotuję zarówno dla okrojonego dataframu, jak i tego z imputacją.

cat(paste('Częstości występowania wartości w zmiennej region:\n'))

## Częstości występowania wartości w zmiennej region:

table(new_df$region)

## 
##          BD/ND/FALA        Dolnośląskie  Kujawsko-Pomorskie           Lubelskie 
##                   0                2596                2068                2620 
##            Lubuskie             Łódzkie         Małopolskie         Mazowieckie 
##                1210                2552                2917                4478 
##            Opolskie        Podkarpackie           Podlaskie           Pomorskie 
##                1277                2535                1512                2201 
##             Śląskie      Świętokrzyskie Warmińsko-Mazurskie       Wielkopolskie 
##                4185                1859                1814                3106 
## Zachodnio-Pomorskie 
##                1531

cat('\n')    

mazowieckie <- new_df_clean$p64_var[new_df_clean$region=='Mazowieckie']
slaskie <- new_df_clean$p64_var[new_df_clean$region=='Śląskie']
wielkopolskie <- new_df_clean$p64_var[new_df_clean$region=='Wielkopolskie']

new_df_clean <- new_df_clean %>%
  filter(region %in% c('Mazowieckie', 'Śląskie', 'Wielkopolskie'))

ggplot(new_df_clean, aes(x=p64_var, color=factor(region, labels=c('Mazowieckie', 'Śląskie', 'Wielkopolskie')), fill=factor(region, labels=c('Mazowieckie', 'Śląskie', 'Wielkopolskie')))) +
  geom_density(alpha=0.2, bw=3000) + 
  labs(title='Porównanie gęstości dochodów względem województw (usunięte nan)', 
       x='Dochód',
       y='Gęstość',
       color='Region',
       fill='Region')+
  theme_minimal()

mazowieckie_density <- density(mazowieckie, bw=3000)
slaskie_density <- density(slaskie, bw=3000)
wielkopolskie_density <- density(wielkopolskie, bw=3000)

mazowieckie_cdf <- approxfun(mazowieckie_density$x, cumsum(mazowieckie_density$y) / sum(mazowieckie_density$y))
slaskie_cdf <- approxfun(slaskie_density$x, cumsum(slaskie_density$y) / sum(slaskie_density$y))
wielkopolskie_cdf <- approxfun(wielkopolskie_density$x, cumsum(wielkopolskie_density$y) / sum(wielkopolskie_density$y))

cat(paste('Jaki procent respondentów z województwa mazowieckiego osiąga dochód netto niższy lub równy 5000: ', mazowieckie_cdf(5000)*100, '%', sep="", '\n'))

## Jaki procent respondentów z województwa mazowieckiego osiąga dochód netto niższy lub równy 5000: 83.3657094206963%

cat(paste('Jaki procent respondentów z województwa śląskiego osiąga dochód netto niższy lub równy 5000: ',slaskie_cdf(5000)*100, '%', sep="",'\n'))

## Jaki procent respondentów z województwa śląskiego osiąga dochód netto niższy lub równy 5000: 83.8719849806908%

cat(paste('Jaki procent respondentów z województwa wielkopolskiego osiąga dochód netto niższy lub równy 5000: ',wielkopolskie_cdf(5000)*100, '%', sep="",'\n'))

## Jaki procent respondentów z województwa wielkopolskiego osiąga dochód netto niższy lub równy 5000: 84.9960248631507%

plot(mazowieckie_density$x, mazowieckie_cdf(mazowieckie_density$x), type = "l", col = "pink",
     xlab = "Value", ylab = "CDF", main = "CDF for Mazowieckie data")

plot(slaskie_density$x, slaskie_cdf(slaskie_density$x), type = "l", col = "lightgreen",
     xlab = "Value", ylab = "CDF", main = "CDF for Śląskie data")

plot(wielkopolskie_density$x, wielkopolskie_cdf(wielkopolskie_density$x), type = "l", col = "lightblue",
     xlab = "Value", ylab = "CDF", main = "CDF for Wielkopolskie data")

mazowieckie <- new_df_impute$p64_var[new_df_impute$region=='Mazowieckie']
slaskie <- new_df_impute$p64_var[new_df_impute$region=='Śląskie']
wielkopolskie <- new_df_impute$p64_var[new_df_impute$region=='Wielkopolskie']

new_df_impute <- new_df_impute %>%
  filter(region %in% c('Mazowieckie', 'Śląskie', 'Wielkopolskie'))

ggplot(new_df_impute, aes(x=p64_var, color=factor(region, labels=c('Mazowieckie', 'Śląskie', 'Wielkopolskie')), fill=factor(region, labels=c('Mazowieckie', 'Śląskie', 'Wielkopolskie')))) +
  geom_density(alpha=0.2, bw=3000) + 
  labs(title='Porównanie gęstości dochodów względem województw (imputacja nan)', 
       x='Dochód',
       y='Gęstość',
       color='Region',
       fill='Region')+
  theme_minimal()

mazowieckie_density_impute <- density(mazowieckie, bw=3000)
slaskie_density_impute <- density(slaskie, bw=3000)
wielkopolskie_density_impute <- density(wielkopolskie, bw=3000)

mazowieckie_cdf_impute <- approxfun(mazowieckie_density_impute$x, cumsum(mazowieckie_density_impute$y) / sum(mazowieckie_density_impute$y))
slaskie_cdf_impute <- approxfun(slaskie_density_impute$x, cumsum(slaskie_density_impute$y) / sum(slaskie_density_impute$y))
wielkopolskie_cdf_impute <- approxfun(wielkopolskie_density_impute$x, cumsum(wielkopolskie_density_impute$y) / sum(wielkopolskie_density_impute$y))

cat(paste('Jaki procent respondentów z województwa mazowieckiego osiąga dochód netto niższy lub równy 5000: ', mazowieckie_cdf_impute(5000)*100, '%', sep="", '\n'))

## Jaki procent respondentów z województwa mazowieckiego osiąga dochód netto niższy lub równy 5000: 84.7378227280673%

cat(paste('Jaki procent respondentów z województwa śląskiego osiąga dochód netto niższy lub równy 5000: ',slaskie_cdf_impute(5000)*100, '%', sep="",'\n'))

## Jaki procent respondentów z województwa śląskiego osiąga dochód netto niższy lub równy 5000: 85.1981556791398%

cat(paste('Jaki procent respondentów z województwa wielkopolskiego osiąga dochód netto niższy lub równy 5000: ',wielkopolskie_cdf_impute(5000)*100, '%', sep="",'\n'))

## Jaki procent respondentów z województwa wielkopolskiego osiąga dochód netto niższy lub równy 5000: 85.693087036592%

plot(mazowieckie_density_impute$x, mazowieckie_cdf_impute(mazowieckie_density_impute$x), type = "l", col = "pink",
     xlab = "Value", ylab = "CDF", main = "CDF for Mazowieckie data")

plot(slaskie_density_impute$x, slaskie_cdf_impute(slaskie_density_impute$x), type = "l", col = "lightgreen",
     xlab = "Value", ylab = "CDF", main = "CDF for Śląskie data")

plot(wielkopolskie_density_impute$x, wielkopolskie_cdf_impute(wielkopolskie_density_impute$x), type = "l", col = "lightblue",
     xlab = "Value", ylab = "CDF", main = "CDF for Wielkopolskie data")

Wnioski

Działając na zbiorze danych Diagnoza, a właściwie na jego części (3 zmienne), chcąc oszacować rozkład dochodu miesięcznego netto według płci i województw, osiągnięty został zamierzony efekt. To znaczy, podejście do pracy miało mieć charakter holistyczny, to znaczy przemyślany i zwracający uwagę na szersze spektrum. Do braków danych podszedłem na dwa sposoby, skutecznie implementując podejście dla każdego z wątków pracy.

Główne wyniki po pracy:

dochód względem płci:

• 99.6% procent respondentów (kobiet) osiąga dochód netto niższy lub równy 10000 (dane po usunięciu NAN),
• 99.4% procent respondentów (kobiet) osiąga dochód netto niższy lub równy 10000 (dane imputowane),
• 99,3% procent respondentów (mężczyzn) osiąga dochód netto niższy lub równy 10000 (dane po usunięciu NAN),
• 98,8% procent respondentów (mężczyzn) osiąga dochód netto niższy lub równy 10000 (dane imputowane).

dochód względem wybranych województw:

• 83.4% procent respondentów z województwa mazowieckiego osiąga dochód netto niższy lub równy 5000 (dane po usunięciu NAN),
• 84.7% procent respondentów z województwa mazowieckiego osiąga dochód netto niższy lub równy 5000 (dane imputowane),
• 83,9% procent respondentów z województwa śląskiego osiąga dochód netto niższy lub równy 5000 (dane po usunięciu NAN),
• 85,2% procent respondentów z województwa śląskiego osiąga dochód netto niższy lub równy 5000 (dane imputowane),
• 85% procent respondentów województwa wielkopolskiego osiąga dochód netto niższy lub równy 5000 (dane po usunięciu NAN),
• 85,7% procent respondentów województwa wielkopolskiego osiąga dochód netto niższy lub równy 5000 (dane imputowane).

Zatem warto zauważyć, że nieznaczna część respondentów osiągała dochód powyżej 10000 złotych miesięcznie na rękę jeżeli chodzi o porównanie względem płci. Uogólniając zaledwie 1 osoba na 100 bez podziału na płeć zarabia powyżej 10000 złotych. Natomiast względem województw obniżenie poziomu o połowę do 5000 złotych nie spowodowało znaczącego spadku w liczności respondentów deklarujących taki poziom dochodów. Można zatem uogólnić, że ponad 8 na 10 osób z wybranych województw osiąga dochód miesięczny netto nieprzekraczający 5000 złotych.