#Загружаем библиотеки:
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.