#Загружаем библиотеки:
library(readr)
library(DT) #для построения таблиц
library(dplyr)
library(stringr)
library(psych) #загрузим библиотеку psych и с ее помощью проанализируем описательные статистики (попутно добавляя те, которые не предусмотрены функцией describe и удаляя лишние)
library(knitr)
library(ggplot2)
library(lubridate)
library(readxl)
library(rio)
library(tidyverse)
library(foreign)
library(coin)
options(scipen=99999)

#Функция для нахождения моды:
mode <- function(x) {
d <- density(x)
d$x[which.max(d$y)]
}
convert("pay_gap.dta", "RLMS_2019.sav")
data = read.spss("RLMS_2019.sav",to.data.frame=T)

В первую очередь посчитаем не скорректированный pay gap. Для этого воспользуемся переменной j13_2 - ответ на вопрос “За последние 12 месяцев какова была Ваша среднемесячная зарплата на этом предприятии (на основной работе) после вычета налогов?”

data_gap = data %>% select(h5,j13_2, j6_2) %>% rename(gender=h5, mwage=j13_2, hours = j6_2) %>% filter(mwage!="ЗАТРУДНЯЮСЬ ОТВЕТИТЬ" & mwage!="ОТКАЗ ОТ ОТВЕТА" & mwage!="НЕТ ОТВЕТА") 
data_gap$mwage=as.numeric(levels(data_gap$mwage))[data_gap$mwage]
c = data_gap %>% group_by(gender) %>% summarise(median = median(mwage))
c
## # A tibble: 2 x 2
##   gender  median
##   <fct>    <dbl>
## 1 МУЖСКОЙ  30000
## 2 ЖЕНСКИЙ  20000
gap1 = 1-c[2,2]/c[1,2]
gap1
##      median
## 1 0.3333333

Gap составил 33%. Проверим статистическим тестом, действительно ли распределение зп зависит от пола: Я проверил по Манна-Уитни (он вроде надежнее t-теста), но, по-моему, все тесты примерно одинаковый результат будут выдавать.

wilcox.test(mwage ~ gender, data = data_gap)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  mwage by gender
## W = 7576378, p-value < 0.00000000000000022
## alternative hypothesis: true location shift is not equal to 0

Зависимость очевидно есть, посмотрим, как она выражается на боксплотах:

ggplot(data_gap %>% filter(mwage<150000))+
  geom_boxplot(aes(x = gender, y = mwage/1000))

Посмотрим на описательные статистики:

female=data_gap %>% filter(gender=="ЖЕНСКИЙ")
male=data_gap %>% filter(gender=="МУЖСКОЙ")
des= describe(male$mwage)  %>% select(-vars, -n, -sd, -trimmed, -mad, -se)%>% rename(asimetria = skew) %>% mutate(firstQu = quantile(male$mwage,0.25)) %>% mutate(thirdQu = quantile(male$mwage,0.75)) %>% mutate(sko = sd(male$mwage)) %>% mutate(moda = mode(male$mwage)) %>% mutate_all(funs(ifelse(.<1, round(.,3), ifelse(.<10,round(.,2),ifelse(.<100,round(.,1),round(.,0)))))) 
kable(des, row.names = F, caption = "Описательные характеристики показателей зарплат у мужчин")
Описательные характеристики показателей зарплат у мужчин
mean median min max range asimetria kurtosis firstQu thirdQu sko moda
34960 30000 0 420000 420000 4.41 47.5 20000 40000 23813 27770
des1= describe(female$mwage)  %>% select(-vars, -n, -sd, -trimmed, -mad, -se)%>% rename(asimetria = skew) %>% mutate(firstQu = quantile(female$mwage,0.25)) %>% mutate(thirdQu = quantile(female$mwage,0.75)) %>% mutate(sko = sd(female$mwage)) %>% mutate(moda = mode(female$mwage)) %>% mutate_all(funs(ifelse(.<1, round(.,3), ifelse(.<10,round(.,2),ifelse(.<100,round(.,1),round(.,0))))))
kable(des1, row.names = F, caption = "Описательные характеристики показателей зарплат у женщин")
Описательные характеристики показателей зарплат у женщин
mean median min max range asimetria kurtosis firstQu thirdQu sko moda
24661 20000 0 180000 180000 2.52 12.1 15000 30000 16394 15526

Сразу видно большую разницу в максимумах и в мерах центральности.

И заодно посмотрим на графики:

ggplot(male) +geom_histogram(aes(x = mwage/1000), fill='pink',col='#483D8B',binwidth = 10)+
  ggtitle("Распределение еженедельных зарплат мужчин")+
  xlab("Зарплата в тыс.рублей")+
  ylab("Абсолютная частота")

ggplot(male)+geom_density(aes(x=mwage/1000))+
  ggtitle("График эмпирической плотности распределения зарплат\nмужчин")+
  ylab("")+
  xlab("Зарплата в тыс.рублей")+
  geom_vline(xintercept = mean(male$mwage/1000), color = 2, linetype = "dashed")+
  annotate("text", x = mean(male$mwage/1000),y=0.024, label = " Среднее", color = 2)+
  geom_vline(xintercept = mode(male$mwage/1000), color = 3, linetype = "dashed")+
  annotate("text", x = mode(male$mwage/1000), y = 0.017, label = "Мода", color = 3)+
  geom_vline(xintercept = median(male$mwage/1000), color = 4, linetype = "dashed")+
  annotate("text", x = median(male$mwage/1000), y = 0.02, label = "Медиана", color = 4)

ggplot(female) +geom_histogram(aes(x = female$mwage/1000),fill = "#87CEEB", col="black", binwidth = 10)+
  ggtitle("Распределение еженедельных зарплат женщин")+
  xlab("Зарплата в тыс.рублей")+
  ylab("Абсолютная частота")

ggplot(female)+geom_density(aes(x=female$mwage/1000))+
  ggtitle("График эмпирической плотности распределения зарплат\nженщин")+
  ylab("")+
  xlab("Зарплата в тыс.рублей")+
  geom_vline(xintercept = mean(female$mwage/1000), color = 2, linetype = "dashed")+
  annotate("text", x = mean(female$mwage/1000),y=0.041, label = " Среднее", color = 2)+
  geom_vline(xintercept = mode(female$mwage/1000), color = 3, linetype = "dashed")+
  annotate("text", x = mode(female$mwage/1000), y = 0.021, label = "Мода", color = 3)+
  geom_vline(xintercept = median(female$mwage/1000), color = 4, linetype = "dashed")+
  annotate("text", x = median(female$mwage/1000), y = 0.029, label = "Медиана", color = 4)

Распределения вцелом схожи, но зп женщин сильно сдвинута (нужно поумнее как-то это написать).

Теперь попробуем хоть как-то скорректировать наш gap. Возьмем не среднемесячную, а почасовую оплату:

data_gap = data_gap %>% filter(hours != "ЗAТРУДНЯЮCЬ ОТВЕТИТЬ" & hours != "ОТКAЗ ОТ ОТВЕТA" & hours != "НЕТ ОТВЕТA")
data_gap$hours=as.numeric(levels(data_gap$hours))[data_gap$hours]
data_gap$hours=data_gap$hours*4
data_gap$hwage = data_gap$mwage/(data_gap$hours)
b = data_gap %>% group_by(gender) %>% summarise(median = median(hwage))
b
## # A tibble: 2 x 2
##   gender  median
##   <fct>    <dbl>
## 1 МУЖСКОЙ   167.
## 2 ЖЕНСКИЙ   125
gap2 = 1-b[2,2]/b[1,2]
gap2
##   median
## 1   0.25

Gap теперь 25% - меньше, видимо значимая доля pay gap появляется из-за того что женщины просто работают меньше. Проверим статтестом зависимость

wilcox.test(hwage ~ gender, data = data_gap)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  hwage by gender
## W = 6294394, p-value < 0.00000000000000022
## alternative hypothesis: true location shift is not equal to 0

Зависит (Нужно ли строить графики?)

Ну чего, посмотрим теперь хоть как-то вертикальную иерархию. Изучим интересную переменную “Есть ли у вас подчиненные”

data_vert = data %>% select(h5,j6, j6_0) %>% rename(gender = h5, subord = j6,num_sub = j6_0) %>% filter(subord != "ЗAТРУДНЯЮCЬ ОТВЕТИТЬ" & subord != "ОТКAЗ ОТ ОТВЕТA" & subord != "НЕТ ОТВЕТA" & subord != "ЗАТРУДНЯЮСЬ ОТВЕТИТЬ" & subord != "ОТКАЗ ОТ ОТВЕТА" & subord != "НЕТ ОТВЕТА")
data_vert$subord %>% summary()
##                   Да                  Нет ЗАТРУДНЯЮСЬ ОТВЕТИТЬ 
##                 1497                 6336                    0 
##      ОТКАЗ ОТ ОТВЕТА           НЕТ ОТВЕТА 
##                    0                    0
ggplot(data_vert)+
  geom_bar(aes(x = gender,fill = subord),color = "black")

Не очень… Ситуация примерно одинаковая. Посмотрим, может, по количеству подчиненных будет интереснее картина

data_vert= data_vert %>% filter(num_sub != "ЗAТРУДНЯЮCЬ ОТВЕТИТЬ" & num_sub != "ОТКAЗ ОТ ОТВЕТA" & num_sub != "НЕТ ОТВЕТA")
data_vert$num_sub=as.numeric(levels(data_vert$num_sub))[data_vert$num_sub]
data_vert$num_sub1 = case_when(
  data_vert$num_sub < 10 ~ "<10",
  data_vert$num_sub < 30 ~ "10 - 30",
  data_vert$num_sub < 60 ~ "30-60",
  data_vert$num_sub < 100 ~ "60-100",
  TRUE ~ ">100"
)
ggplot(data_vert)+
  geom_bar(aes(x = gender,fill = num_sub1),color = "black")

Большая часть - это те, у которых меньше 10 подчиненных. Уберем их

data_vert = data_vert %>% filter(num_sub>10)
data_vert$num_sub2 = case_when(
  data_vert$num_sub < 25 ~ "<25",
  data_vert$num_sub < 60 ~ "25 - 60",
  data_vert$num_sub < 125 ~ "60-125",
  data_vert$num_sub < 250 ~ "125-250",
  data_vert$num_sub < 500 ~ "250-500",
  TRUE ~ ">500"
)
ggplot(data_vert)+
  geom_bar(aes(x = gender,fill = num_sub2),color = "black")

Важно в первую очередь то, что мужчин в этой группе в принципе больше. Ну и в каждой отдельной группе видно их преимущество (кроме 250-500…)

Давайте смотреть профессии. Наша гипотеза в том, что wage gap объясняется тем, что женщин банально меньше на высокооплачиваемых должностях.

data_occ = data %>% select(occup08,h5,j13_2) %>% rename(gender = h5, occ = occup08,mwage = j13_2) %>% filter(!is.na(occ)) %>% filter(occ != "ЗАТРУДНЯЮСЬ ОТВЕТИТЬ" & occ != "ОТКАЗ ОТ ОТВЕТА" & occ != "НЕТ ОТВЕТА")
ggplot(data_occ)+
  geom_bar(aes(x=occ, fill = gender))+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

На графике обратная картина……

ggplot(data_occ)+
  geom_bar(aes(x=gender, fill = occ),position = "fill")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Можем обвинить данные, потому что не все дают ответ на вопрос про род деятельности. Ну и +наверное деление на такие группы не совсем точное

Посмотрим на pay gap в каждом из видов деятельности:

data_occ = data_occ %>% filter(!is.na(mwage)) %>% filter(mwage != "ЗАТРУДНЯЮСЬ ОТВЕТИТЬ" & mwage != "ОТКАЗ ОТ ОТВЕТА" & mwage != "НЕТ ОТВЕТА")
data_occ$mwage=data_occ$mwage=as.numeric(levels(data_occ$mwage))[data_occ$mwage]
d = data_occ %>% group_by(occ,gender) %>% summarise(median = median(mwage))
d
## # A tibble: 20 x 3
## # Groups:   occ [10]
##    occ                                                             gender median
##    <fct>                                                           <fct>   <dbl>
##  1 военнослужащие                                                  МУЖСК…  40900
##  2 военнослужащие                                                  ЖЕНСК…  40000
##  3 законодатели; крупные чиновники; руководители высш. и сред. зв… МУЖСК…  40000
##  4 законодатели; крупные чиновники; руководители высш. и сред. зв… ЖЕНСК…  30000
##  5 специалисты высшего уровня квалификации                         МУЖСК…  40000
##  6 специалисты высшего уровня квалификации                         ЖЕНСК…  25000
##  7 специалисты среднего уровня квалификации; чиновники             МУЖСК…  35000
##  8 специалисты среднего уровня квалификации; чиновники             ЖЕНСК…  22000
##  9 служащие офисные и по обслуживанию клиентов                     МУЖСК…  32000
## 10 служащие офисные и по обслуживанию клиентов                     ЖЕНСК…  20000
## 11 работники сферы торговли и услуг                                МУЖСК…  22000
## 12 работники сферы торговли и услуг                                ЖЕНСК…  17000
## 13 квалифицированные работники сельского, лесного хоз-ва и рыбово… МУЖСК…  27000
## 14 квалифицированные работники сельского, лесного хоз-ва и рыбово… ЖЕНСК…  47000
## 15 квалифицированные рабочие, занятые ручным трудом                МУЖСК…  30000
## 16 квалифицированные рабочие, занятые ручным трудом                ЖЕНСК…  22000
## 17 квалифицированные рабочие, использующие машины и механизмы      МУЖСК…  30000
## 18 квалифицированные рабочие, использующие машины и механизмы      ЖЕНСК…  20000
## 19 неквалифицированные рабочие всех отраслей                       МУЖСК…  20000
## 20 неквалифицированные рабочие всех отраслей                       ЖЕНСК…  14000
ggplot(d)+
  geom_bar(aes(x = occ,y = median, fill = gender), position = "dodge", stat = "identity")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  coord_flip()

Почти во всех группах наблюдается pay gap.