проект группы 21
2022-12-16
#Код 21 группы к проекту по эконометрике
Подготовка к работе
Загрузка пакетов
library(readr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(GGally)
library(sandwich)
library(lmtest)
library(car)
library(broom)
library(xtable)
library(ggpubr)
library(stargazer)
library(modelsummary)
library(haven)
library(readr)
library(margins)Установка рабочей директории
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))Собственные функции
number <- function(x, na.rm = TRUE){return(sum(!is.na(x)))}Загрузка данных
data_0 <- read_sav("r30iall_53.sav")
data_0## # A tibble: 17,605 × 695
## idind zredid_i zid_i zid_h yid_i yid_h xid_i xid_h wid_i wid_h vid_i
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 55501 1000101 10001 1000101 10001 1.00e6 10001 1000101 10001 1000101
## 2 3 56001 1001201 10012 1001201 10012 1.00e6 10012 1001201 10012 1001201
## 3 6 70204 1501404 15014 1501404 15014 1.50e6 15014 1501404 15014 1501404
## 4 7 55601 1000301 10003 1000301 10003 1.00e6 10003 1000301 10003 1000301
## 5 9 55701 1000401 10004 1000401 10004 1.00e6 10004 NA NA 1000401
## 6 28 56501 1002101 10021 1002101 10021 1.00e6 10021 1002101 10021 NA
## 7 29 70502 1502002 15020 1502002 15020 1.50e6 15020 1502002 15020 1002102
## 8 30 56503 1002103 10021 1002103 10021 1.00e6 10021 1002103 10021 1002103
## 9 36 56701 1002401 10024 1002401 10024 1.00e6 10024 1002401 10024 NA
## 10 37 70402 1501902 15019 1501902 15019 1.50e6 15019 1501902 15019 1002402
## # … with 17,595 more rows, and 684 more variables: vid_h <dbl>, uid_i <dbl>,
## # uid_h <dbl>, tid_i <dbl>, tid_h <dbl>, sid_i <dbl>, sid_h <dbl>,
## # rid_i <dbl>, rid_h <dbl>, qid_i <dbl>, qid_h <dbl>, pid_i <dbl>,
## # pid_h <dbl>, oid_i <dbl>, oid_h <dbl>, nid_i <dbl>, nid_h <dbl>,
## # mid_i <dbl>, mid_h <dbl>, lid_i <dbl>, lid_h <dbl>, kid_i <dbl>,
## # kid_h <dbl>, jid_i <dbl>, jid_h <dbl>, iid_i <dbl>, iid_h <dbl>,
## # hid_i <dbl>, hid_h <dbl>, gid_i <dbl>, gid_h <dbl>, fid_i <dbl>, …Чистка данных
Значение переменных в базе данных:
status - тип населенного пункта
z_marst - СЕМЕЙНОЕ ПОЛОЖЕНИЕ В ТЕКУЩЕЙ ВОЛНЕ
z_occup08 - ПРОФЕССИОНАЛЬНАЯ ГРУППА
z_diplom - ЗАКОНЧЕННОЕ ОБРАЗО (ГРУППА)
z_age - Количество полных лет (возраст респондента)
zh5 - Пол респондента
zj13.2 - За последние 12 месяцев какова была Ваша среднемесячная зарплата на этом предприятии после вычета налогов - независимо от того, платят Вам ее вовремя или нет?
zj161.3y - Сколько полных лет и месяцев составляет Ваш официально оформленный общий трудовой 99999996 стаж, не считая времени учебы на дневном отделении ВУЗа или техникума? (лет)
zj90 - Какой ответ лучше всего описывает Ваше основное занятие в настоящее время?
data_0 = data_0 %>% select(status,z_marst,z_occup08,z_diplom,z_age,zh5,
zj13.2,zj161.3y,zj90) %>% na.omit()
data_0 = data_0 %>% filter(status < 99999990,z_marst < 99999990,
zj13.2< 99999990,z_age< 99999990,z_occup08<99999990,
z_diplom<99999990,zj161.3y< 99999990,zh5<99999990,
zj90<14,zj90>9)
#data=data %>% select(-zj90)Оставляет только ЭАН и фильтруем по роду деятельности
data = data_0 %>% filter(z_age >= 16 & z_age <= 72) %>% filter(zj90 %in% c(10, 11, 12, 13))Название колонок
colnames(data) <- c("place","family","proff","educ","age","sex","wage",
"work_exper","Пол", "образование", "место",
"семейное полож", "Работа", "working_hours")Cмещенность выборки
Результат t-теста на разницу в средних значени зп до фильтрации выборки по возрасту и по виду занятости. P-value > 0,05, значит нулевая гипотеза об истинном равенстве средних не отвергается. И выборка не смещенная
t.test(data$wage, data_0$zj13.2)##
## Welch Two Sample t-test
##
## data: data$wage and data_0$zj13.2
## t = 0.1068, df = 11762, p-value = 0.915
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -782.9196 873.1496
## sample estimates:
## mean of x mean of y
## 34266.49 34221.38Разведывательный анализ
График распределения заработных плат индивидов
ggplot(data, aes(wage)) +
geom_histogram(aes(y = ..density..), color = "black", alpha = 0.25) +
geom_density(color = "blue", adjust = 4)+
ggtitle("Распределение заработных плат индивидов")+
xlab('Среднемесячная зарплата')+
ylab('Плотность')+
scale_x_continuous(breaks = seq(0, 400000, by = 50000))## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Логарифмирование заработной платы и фильтрация нулевых значений
data$log_wage<- log(data$wage)
data = data %>% filter(log_wage > 0)График распределения прологарифмированных зарплат
ggplot(data, aes(log_wage)) +
geom_histogram(aes(y = ..density..), color = "black", alpha = 0.25) +
geom_density(color = "blue", adjust = 4)+
ggtitle("Распределение заработных плат (log) индивидов")+
xlab('Среднемесячная зарплата')+
ylab('Плотность')+
scale_x_continuous(breaks = seq(0, 13, by = 1))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Бинарная переменная big_city
data$big_city = ifelse(data$place %in% c(1, 2), 1, 0)
data$big_city =factor(data$big_city,c(0,1), labels=c("Маленький н.п.", "Большой н.п."))
data$place =as.factor(data$place)
data$place =factor(data$place,c(1,2,3,4), labels=c("Обл. центр", "Город","ПГТ", "Село"))График распределения по типу н.п.
ggplot(data, aes(x = big_city, fill = place)) +
geom_histogram(stat = 'count', position = 'dodge') +
ggtitle("Распределение респондентов по типу населенного пункта")+
xlab('Р населённого пункта')+
ylab('Количество')+
scale_fill_discrete(name = "Тип н.п.")## Warning in geom_histogram(stat = "count", position = "dodge"): Ignoring unknown
## parameters: `binwidth`, `bins`, and `pad`
По числам
data %>% group_by(big_city, place) %>% summarize(count = n())## `summarise()` has grouped output by 'big_city'. You can override using the
## `.groups` argument.## # A tibble: 4 × 3
## # Groups: big_city [2]
## big_city place count
## <fct> <fct> <int>
## 1 Маленький н.п. ПГТ 354
## 2 Маленький н.п. Село 1166
## 3 Большой н.п. Обл. центр 2618
## 4 Большой н.п. Город 1710summary(data$age)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 18.00 34.00 42.00 42.85 51.00 72.00data$educ =as.factor(data$educ)
data$sex =as.factor(data$sex)
data$family =as.factor(data$family)
data$proff =as.factor(data$proff)
data$educ =factor(data$educ,c(1,2,3,4,5,6),labels=c("Окончил 0-6 кл", "Незаконченное сред.обр(7-8)","Незаконченное .обр(7-8+)", "среднее образование", "Среднее спец.образование", "Высшее образование"))
data$sex =factor(data$sex,c(1,2), labels=c("Мужчина","Женщина"))
data$proff =factor(data$proff,c(0,1,2,3,4,5,6,7,8,9), labels=c("военнослужащие","Руководитель", "Высшее специалисты", "Среднее специалисты", "Обслуживание","Услуги","Квалифицированные","Квалифицированные2","Квалифицированные3","Неквалифицированные"))
data$family= factor(data$family,c(1,2,3,4,5,6), labels=c("Н в браке", "Брак", "Живёте вместе", "Разведены","Вдова(вдовец)", "Оф.зарегистрированы"))График распределения по полу
ggplot(data, aes(x = big_city, fill = sex)) +
geom_histogram(stat = 'count', position = 'dodge') +
ggtitle("Распределение респондентов по полу и размеру н. п.")+
xlab('Размер населённого пункта')+
ylab('Количество')+
scale_fill_discrete(name = "Пол")## Warning in geom_histogram(stat = "count", position = "dodge"): Ignoring unknown
## parameters: `binwidth`, `bins`, and `pad`
Описательные статистики
data %>% group_by(big_city, sex) %>% summarize(count = n())## `summarise()` has grouped output by 'big_city'. You can override using the
## `.groups` argument.## # A tibble: 4 × 3
## # Groups: big_city [2]
## big_city sex count
## <fct> <fct> <int>
## 1 Маленький н.п. Мужчина 702
## 2 Маленький н.п. Женщина 818
## 3 Большой н.п. Мужчина 2025
## 4 Большой н.п. Женщина 2303data %>% group_by(big_city, educ) %>% summarize(count = n())## `summarise()` has grouped output by 'big_city'. You can override using the
## `.groups` argument.## # A tibble: 12 × 3
## # Groups: big_city [2]
## big_city educ count
## <fct> <fct> <int>
## 1 Маленький н.п. Окончил 0-6 кл 1
## 2 Маленький н.п. Незаконченное сред.обр(7-8) 33
## 3 Маленький н.п. Незаконченное .обр(7-8+) 130
## 4 Маленький н.п. среднее образование 493
## 5 Маленький н.п. Среднее спец.образование 448
## 6 Маленький н.п. Высшее образование 415
## 7 Большой н.п. Окончил 0-6 кл 3
## 8 Большой н.п. Незаконченное сред.обр(7-8) 38
## 9 Большой н.п. Незаконченное .обр(7-8+) 206
## 10 Большой н.п. среднее образование 1115
## 11 Большой н.п. Среднее спец.образование 1142
## 12 Большой н.п. Высшее образование 1824data %>% group_by(big_city, family) %>% summarize(count = n())## `summarise()` has grouped output by 'big_city'. You can override using the
## `.groups` argument.## # A tibble: 12 × 3
## # Groups: big_city [2]
## big_city family count
## <fct> <fct> <int>
## 1 Маленький н.п. Н в браке 206
## 2 Маленький н.п. Брак 919
## 3 Маленький н.п. Живёте вместе 180
## 4 Маленький н.п. Разведены 123
## 5 Маленький н.п. Вдова(вдовец) 85
## 6 Маленький н.п. Оф.зарегистрированы 7
## 7 Большой н.п. Н в браке 627
## 8 Большой н.п. Брак 2507
## 9 Большой н.п. Живёте вместе 544
## 10 Большой н.п. Разведены 436
## 11 Большой н.п. Вдова(вдовец) 191
## 12 Большой н.п. Оф.зарегистрированы 23summary(data$work_exper)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 10.00 19.00 19.78 28.00 57.00data %>% group_by(big_city, proff) %>% summarize(count = n())## `summarise()` has grouped output by 'big_city'. You can override using the
## `.groups` argument.## # A tibble: 20 × 3
## # Groups: big_city [2]
## big_city proff count
## <fct> <fct> <int>
## 1 Маленький н.п. военнослужащие 7
## 2 Маленький н.п. Руководитель 88
## 3 Маленький н.п. Высшее специалисты 266
## 4 Маленький н.п. Среднее специалисты 299
## 5 Маленький н.п. Обслуживание 76
## 6 Маленький н.п. Услуги 261
## 7 Маленький н.п. Квалифицированные 5
## 8 Маленький н.п. Квалифицированные2 160
## 9 Маленький н.п. Квалифицированные3 230
## 10 Маленький н.п. Неквалифицированные 128
## 11 Большой н.п. военнослужащие 22
## 12 Большой н.п. Руководитель 294
## 13 Большой н.п. Высшее специалисты 855
## 14 Большой н.п. Среднее специалисты 945
## 15 Большой н.п. Обслуживание 272
## 16 Большой н.п. Услуги 618
## 17 Большой н.п. Квалифицированные 3
## 18 Большой н.п. Квалифицированные2 559
## 19 Большой н.п. Квалифицированные3 500
## 20 Большой н.п. Неквалифицированные 260Регрессионный анализ
Модель 1
reg1 = lm(data = data, log_wage ~ big_city+sex+family+age+educ+proff)
cov1 <- vcovHC(reg1, type = "HC0")
coeftest(reg1, df = Inf, vcov = vcovHC, type = "HC0")##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 10.48744200 0.19200610 54.6204 < 2.2e-16 ***
## big_cityБольшой н.п. 0.16956733 0.01492658 11.3601 < 2.2e-16 ***
## sexЖенщина -0.32573085 0.01501165 -21.6985 < 2.2e-16 ***
## familyБрак 0.08887664 0.02066689 4.3004 1.705e-05 ***
## familyЖивёте вместе 0.09366830 0.02540738 3.6867 0.0002272 ***
## familyРазведены 0.12676210 0.02882583 4.3975 1.095e-05 ***
## familyВдова(вдовец) 0.11374537 0.03730841 3.0488 0.0022977 **
## familyОф.зарегистрированы 0.11356748 0.07014938 1.6189 0.1054607
## age -0.00744394 0.00065813 -11.3107 < 2.2e-16 ***
## educНезаконченное сред.обр(7-8) 0.06072795 0.19041918 0.3189 0.7497893
## educНезаконченное .обр(7-8+) 0.12790805 0.18247526 0.7010 0.4833273
## educсреднее образование 0.14465956 0.18103272 0.7991 0.4242442
## educСреднее спец.образование 0.15020258 0.18130618 0.8284 0.4074173
## educВысшее образование 0.33854033 0.18151076 1.8651 0.0621638 .
## proffРуководитель 0.16359243 0.06383581 2.5627 0.0103860 *
## proffВысшее специалисты -0.02774956 0.06023803 -0.4607 0.6450388
## proffСреднее специалисты -0.09282023 0.05950758 -1.5598 0.1188059
## proffОбслуживание -0.16174949 0.06314988 -2.5614 0.0104264 *
## proffУслуги -0.32587791 0.06033106 -5.4015 6.609e-08 ***
## proffКвалифицированные -0.05055512 0.21189717 -0.2386 0.8114287
## proffКвалифицированные2 -0.17617646 0.06002036 -2.9353 0.0033325 **
## proffКвалифицированные3 -0.10038233 0.06041901 -1.6614 0.0966259 .
## proffНеквалифицированные -0.46348085 0.06307236 -7.3484 2.006e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Модель 2
reg2 = lm(data = data, log_wage ~ big_city+sex+family+age+educ+proff+big_city*sex)
cov2 <- vcovHC(reg2,type ="HC0", data= data)
coeftest(reg2, df = Inf, vcov = vcovHC, type = "HC0")##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 10.47382811 0.19405732 53.9729 < 2.2e-16 ***
## big_cityБольшой н.п. 0.18370280 0.02198939 8.3542 < 2.2e-16 ***
## sexЖенщина -0.30597123 0.02739342 -11.1695 < 2.2e-16 ***
## familyБрак 0.08849928 0.02067903 4.2797 1.872e-05 ***
## familyЖивёте вместе 0.09335720 0.02541939 3.6727 0.000240 ***
## familyРазведены 0.12639945 0.02883764 4.3831 1.170e-05 ***
## familyВдова(вдовец) 0.11300639 0.03733379 3.0269 0.002471 **
## familyОф.зарегистрированы 0.11528391 0.07034350 1.6389 0.101240
## age -0.00743770 0.00065784 -11.3063 < 2.2e-16 ***
## educНезаконченное сред.обр(7-8) 0.06533813 0.19184977 0.3406 0.733428
## educНезаконченное .обр(7-8+) 0.13135728 0.18389130 0.7143 0.475029
## educсреднее образование 0.14809275 0.18246841 0.8116 0.417017
## educСреднее спец.образование 0.15282869 0.18269587 0.8365 0.402863
## educВысшее образование 0.34193286 0.18291905 1.8693 0.061579 .
## proffРуководитель 0.16313005 0.06376781 2.5582 0.010522 *
## proffВысшее специалисты -0.02849450 0.06015962 -0.4736 0.635751
## proffСреднее специалисты -0.09305592 0.05942994 -1.5658 0.117393
## proffОбслуживание -0.16155886 0.06309658 -2.5605 0.010452 *
## proffУслуги -0.32592945 0.06025814 -5.4089 6.342e-08 ***
## proffКвалифицированные -0.04942076 0.21085396 -0.2344 0.814687
## proffКвалифицированные2 -0.17651089 0.05994705 -2.9444 0.003235 **
## proffКвалифицированные3 -0.09955036 0.06035159 -1.6495 0.099044 .
## proffНеквалифицированные -0.46423472 0.06298501 -7.3706 1.699e-13 ***
## big_cityБольшой н.п.:sexЖенщина -0.02627050 0.02985026 -0.8801 0.378818
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Результаты оценивания регрессии
modelsummary(models = list(reg1, reg2),
vcov = list(cov1, cov2),
statistic = "std.error",
stars = TRUE,
#gof_omit = ".*",
notes = list("В скобках даны робастные стандартные ошибки"))| Model 1 | Model 2 | |
|---|---|---|
| (Intercept) | 10.487*** | 10.474*** |
| (0.192) | (0.194) | |
| big_cityБольшой н.п. | 0.170*** | 0.184*** |
| (0.015) | (0.022) | |
| sexЖенщина | −0.326*** | −0.306*** |
| (0.015) | (0.027) | |
| familyБрак | 0.089*** | 0.088*** |
| (0.021) | (0.021) | |
| familyЖивёте вместе | 0.094*** | 0.093*** |
| (0.025) | (0.025) | |
| familyРазведены | 0.127*** | 0.126*** |
| (0.029) | (0.029) | |
| familyВдова(вдовец) | 0.114** | 0.113** |
| (0.037) | (0.037) | |
| familyОф.зарегистрированы | 0.114 | 0.115 |
| (0.070) | (0.070) | |
| age | −0.007*** | −0.007*** |
| (0.0007) | (0.0007) | |
| educНезаконченное сред.обр(7-8) | 0.061 | 0.065 |
| (0.190) | (0.192) | |
| educНезаконченное .обр(7-8+) | 0.128 | 0.131 |
| (0.182) | (0.184) | |
| educсреднее образование | 0.145 | 0.148 |
| (0.181) | (0.182) | |
| educСреднее спец.образование | 0.150 | 0.153 |
| (0.181) | (0.183) | |
| educВысшее образование | 0.339+ | 0.342+ |
| (0.182) | (0.183) | |
| proffРуководитель | 0.164* | 0.163* |
| (0.064) | (0.064) | |
| proffВысшее специалисты | −0.028 | −0.028 |
| (0.060) | (0.060) | |
| proffСреднее специалисты | −0.093 | −0.093 |
| (0.060) | (0.059) | |
| proffОбслуживание | −0.162* | −0.162* |
| (0.063) | (0.063) | |
| proffУслуги | −0.326*** | −0.326*** |
| (0.060) | (0.060) | |
| proffКвалифицированные | −0.051 | −0.049 |
| (0.212) | (0.211) | |
| proffКвалифицированные2 | −0.176** | −0.177** |
| (0.060) | (0.060) | |
| proffКвалифицированные3 | −0.100+ | −0.100+ |
| (0.060) | (0.060) | |
| proffНеквалифицированные | −0.463*** | −0.464*** |
| (0.063) | (0.063) | |
| big_cityБольшой н.п. × sexЖенщина | −0.026 | |
| (0.030) | ||
| Num.Obs. | 5848 | 5848 |
| R2 | 0.257 | 0.257 |
| R2 Adj. | 0.254 | 0.254 |
| AIC | 8223.9 | 8225.1 |
| BIC | 8384.1 | 8391.9 |
| Log.Lik. | −4087.942 | −4087.539 |
| RMSE | 0.49 | 0.49 |
| Std.Errors | Custom | Custom |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||
| В скобках даны робастные стандартные ошибки |
Корреляция
cor(data$age, data$work_exper)## [1] 0.9063416Предельные эффекты
margins_summary(reg2)## factor AME SE z p lower upper
## age -0.0074 0.0006 -11.7691 0.0000 -0.0087 -0.0062
## big_cityБольшой н.п. 0.1697 0.0148 11.4820 0.0000 0.1407 0.1986
## educВысшее образование 0.3419 0.2448 1.3966 0.1625 -0.1379 0.8218
## educНезаконченное .обр(7-8+) 0.1314 0.2457 0.5347 0.5928 -0.3501 0.6128
## educНезаконченное сред.обр(7-8) 0.0653 0.2510 0.2603 0.7946 -0.4265 0.5572
## educсреднее образование 0.1481 0.2445 0.6056 0.5448 -0.3312 0.6274
## educСреднее спец.образование 0.1528 0.2447 0.6246 0.5323 -0.3268 0.6324
## familyБрак 0.0885 0.0203 4.3639 0.0000 0.0488 0.1282
## familyВдова(вдовец) 0.1130 0.0374 3.0225 0.0025 0.0397 0.1863
## familyЖивёте вместе 0.0934 0.0258 3.6242 0.0003 0.0429 0.1438
## familyОф.зарегистрированы 0.1153 0.0911 1.2654 0.2057 -0.0633 0.2938
## familyРазведены 0.1264 0.0285 4.4389 0.0000 0.0706 0.1822
## proffВысшее специалисты -0.0285 0.0924 -0.3083 0.7579 -0.2097 0.1527
## proffКвалифицированные -0.0494 0.1951 -0.2533 0.8000 -0.4318 0.3329
## proffКвалифицированные2 -0.1765 0.0929 -1.9007 0.0573 -0.3585 0.0055
## proffКвалифицированные3 -0.0996 0.0929 -1.0714 0.2840 -0.2817 0.0826
## proffНеквалифицированные -0.4642 0.0949 -4.8903 0.0000 -0.6503 -0.2782
## proffОбслуживание -0.1616 0.0949 -1.7022 0.0887 -0.3476 0.0245
## proffРуководитель 0.1631 0.0944 1.7286 0.0839 -0.0218 0.3481
## proffСреднее специалисты -0.0931 0.0921 -1.0102 0.3124 -0.2736 0.0875
## proffУслуги -0.3259 0.0927 -3.5144 0.0004 -0.5077 -0.1442
## sexЖенщина -0.3254 0.0151 -21.6099 0.0000 -0.3549 -0.2959