library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(haven)
## Warning: package 'haven' was built under R version 4.0.4
library(writexl)
## Warning: package 'writexl' was built under R version 4.0.5
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.0.5
library(GGally)
## Warning: package 'GGally' was built under R version 4.0.5
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
data2019 = read_sav("r28iall_42.sav")
data = data2019 %>% select(xl90, xj60, xj13.2, x_age, xj21.3, status, x_age, xh5, xj4.1, xj6.2, xj11.1, xj29, xj72.171, xj72.173, xm20.61, xm20.62, xm20.63, xm20.64, xm20.65, xm20.66, xm20.620, xm20.69, xm20.610, xm20.611, xm20.612, xm20.613, xm20.614, xm20.615, xm20.616, xm20.617, xm20.618, xm20.619, xm20.67, xm20.7, xm39)
data = data %>% filter(is.na(xl90) == FALSE & is.na(xj60) == FALSE & is.na(xj13.2) == FALSE)
data$xl90 = as.numeric(data$xl90)
data$xj60 = as.numeric(data$xj60)
data$xj13.2 = as.numeric(data$xj13.2)
data$x_age = as.numeric(data$x_age)
data$xj21.3 = as.numeric(data$xj21.3)
data$status = as.numeric(data$status)
data$xh5 = as.numeric(data$xh5)
data$xj4.1 = as.numeric(data$xj4.1)
data$xj6.2 = as.numeric(data$xj6.2)
data$xj11.1 = as.numeric(data$xj11.1)
data$xj29 = as.numeric(data$xj29)
data$xj72.171 = as.numeric(data$xj72.171)
data$xj72.173 = as.numeric(data$xj72.173)
data$xm20.61 = as.numeric(data$xm20.61)
data$xm20.62 = as.numeric(data$xm20.62)
data$xm20.63 = as.numeric(data$xm20.63)
data$xm20.64 = as.numeric(data$xm20.64)
data$xm20.65 = as.numeric(data$xm20.65)
data$xm20.66 = as.numeric(data$xm20.66)
data$xm20.620 = as.numeric(data$xm20.620)
data$xm20.69 = as.numeric(data$xm20.69)
data$xm20.610 = as.numeric(data$xm20.610)
data$xm20.611 = as.numeric(data$xm20.611)
data$xm20.612 = as.numeric(data$xm20.612)
data$xm20.613 = as.numeric(data$xm20.613)
data$xm20.614 = as.numeric(data$xm20.614)
data$xm20.615 = as.numeric(data$xm20.615)
data$xm20.616 = as.numeric(data$xm20.616)
data$xm20.617 = as.numeric(data$xm20.617)
data$xm20.618 = as.numeric(data$xm20.618)
data$xm20.619 = as.numeric(data$xm20.619)
data$xm20.67 = as.numeric(data$xm20.67)
data$xm20.7 = as.numeric(data$xm20.7)
data$xm39 = as.numeric(data$xm39)
data[is.na(data$xj72.173), "xj72.173"] <- 0 #про детей моложе 18
data[is.na(data$xm20.614), "xm20.614"] <- 2 #про гинекологические заболевания
data = data %>% filter(xm39 < 9999990) %>% filter(xm20.7 < 9999990) %>% filter(xm20.618 < 9999990) %>% filter(xm20.617 < 9999990) %>% filter(xm20.616 < 9999990) %>% filter(xm20.615 < 9999990) %>% filter(xm20.613 < 9999990) %>% filter(xm20.612 < 9999990) %>% filter(xm20.610 < 9999990) %>% filter(xm20.620 < 9999990) %>% filter(xm20.64 < 9999990) %>% filter(xm20.63 < 9999990) %>% filter(xm20.62 < 9999990) %>% filter(xj72.173 < 9999990) %>% filter(xj72.171 < 9999990) %>% filter(xj4.1 < 9999990) %>% filter(xh5 < 9999990) %>% filter(status < 9999990) %>% filter(x_age < 9999990) %>% filter(xl90 < 9999990) %>% filter(xj13.2 < 9999990) %>% filter(xj60 < 9999990) %>% filter(xj21.3 < 9999990) %>% filter(xj6.2 < 99999990) %>% filter(xj11.1 < 9999990) %>% filter(xj29 < 9999990) %>% filter(xm20.61 < 9999990) %>% filter(xm20.65 < 9999990) %>% filter(xm20.66 < 9999990) %>% filter(xm20.69 < 9999990) %>% filter(xm20.611 < 9999990) %>% filter(xm20.619 < 9999990) %>% filter(xm20.67 < 9999990) %>% filter(xm20.614 < 9999990)
data[data$xm20.61 == 2, "xm20.61"] <- 0
data[data$xm20.62 == 2, "xm20.62"] <- 0
data[data$xm20.63 == 2, "xm20.63"] <- 0
data[data$xm20.64 == 2, "xm20.64"] <- 0
data[data$xm20.65 == 2, "xm20.65"] <- 0
data[data$xm20.66 == 2, "xm20.66"] <- 0
data[data$xm20.67 == 2, "xm20.67"] <- 0
data[data$xm20.69 == 2, "xm20.69"] <- 0
data[data$xm20.610 == 2, "xm20.610"] <- 0
data[data$xm20.611 == 2, "xm20.611"] <- 0
data[data$xm20.612 == 2, "xm20.612"] <- 0
data[data$xm20.613 == 2, "xm20.613"] <- 0
data[data$xm20.614 == 2, "xm20.614"] <- 0
data[data$xm20.615 == 2, "xm20.615"] <- 0
data[data$xm20.616 == 2, "xm20.616"] <- 0
data[data$xm20.617 == 2, "xm20.617"] <- 0
data[data$xm20.618 == 2, "xm20.618"] <- 0
data[data$xm20.619 == 2, "xm20.619"] <- 0
data[data$xm20.620 == 2, "xm20.620"] <- 0
data = data %>% mutate(quant_chron_zabol = xm20.61 + xm20.62 + xm20.620 + xm20.63 + xm20.64 + xm20.65 + xm20.66 + xm20.67 + xm20.69 + xm20.610 + xm20.611 + xm20.612 + xm20.613 + xm20.614 + xm20.615 + xm20.616 + xm20.617 + xm20.618 + xm20.619) %>% dplyr::select(-xm20.61,-xm20.62,-xm20.620,-xm20.63,-xm20.64,-xm20.65,-xm20.66,-xm20.67,-xm20.69,-xm20.610,-xm20.611,-xm20.612,-xm20.613,-xm20.614,-xm20.615,-xm20.616,-xm20.617,-xm20.618,-xm20.619)
data$xj21.3 = as.factor(data$xj21.3)
data$status = as.factor(data$status)
data$xh5 = as.factor(data$xh5)
data$xj4.1 = as.factor(data$xj4.1)
data$xj11.1 = as.factor(data$xj11.1)
data$xj29 = as.factor(data$xj29)
data$xj72.171 = as.factor(data$xj72.171)
data$xj72.173 = as.factor(data$xj72.173)
data$xm20.7 = as.factor(data$xm20.7)
data$xm39 = as.factor(data$xm39)
data = data %>% filter(xj60 > 0, xj13.2 > 0, xj6.2 < 100)
data$xl90 = data$xl90/12
colnames(data)
## [1] "xl90" "xj60" "xj13.2"
## [4] "x_age" "xj21.3" "status"
## [7] "xh5" "xj4.1" "xj6.2"
## [10] "xj11.1" "xj29" "xj72.171"
## [13] "xj72.173" "xm20.7" "xm39"
## [16] "quant_chron_zabol"
colnames(data) = c("sick_days", "monthly_income", "average_wage", "age", "harmful_work", "status", "gender", "branch", "hours", "official_work", "entrepreneurship", "children", "children_count", "disability", "operations", "quant_chron_zabol")
summary(data)
## sick_days monthly_income average_wage age
## Min. : 0.08333 Min. : 1800 Min. : 3000 Min. :18.00
## 1st Qu.: 0.58333 1st Qu.: 20000 1st Qu.: 18000 1st Qu.:33.00
## Median : 0.83333 Median : 30000 Median : 25000 Median :41.00
## Mean : 1.34806 Mean : 40875 Mean : 30641 Mean :42.11
## 3rd Qu.: 1.50000 3rd Qu.: 43000 3rd Qu.: 38000 3rd Qu.:51.00
## Max. :22.91667 Max. :2550000 Max. :420000 Max. :80.00
##
## harmful_work status gender branch hours official_work
## 1: 211 1:639 1:529 14 :242 Min. : 4.00 1:1259
## 2:1096 2:388 2:778 10 :159 1st Qu.:40.00 2: 48
## 3: 57 7 :132 Median :40.00
## 4:223 12 :102 Mean :42.01
## 1 : 87 3rd Qu.:45.00
## 6 : 74 Max. :96.00
## (Other):511
## entrepreneurship children children_count disability operations
## 1: 44 1:1050 0:715 1: 43 1: 111
## 2:1263 2: 257 1:330 2:1264 2:1196
## 2:218
## 3: 39
## 4: 3
## 5: 1
## 6: 1
## quant_chron_zabol
## Min. : 0.000
## 1st Qu.: 0.000
## Median : 1.000
## Mean : 1.744
## 3rd Qu.: 3.000
## Max. :13.000
##
ggplot(data, mapping = aes(x = as.numeric(monthly_income), y = sick_days)) +
geom_point(alpha = 0.5) +
scale_x_log10() +
scale_y_log10() +
ggtitle('Соотношение логарифма количества пропущенных дней на логарифм дохода') +
xlab('Логарифм общего дохода за последний месяц') +
ylab('Логарифм среднего количества пропущенных дней за месяц')
ggplot(data, mapping = aes(x = average_wage , y = sick_days)) +
geom_point(alpha = 0.5) +
scale_x_log10() +
scale_y_log10() +
ggtitle('Соотношение логарифма количества пропущенных дней на логарифм зп') +
xlab('Логарифм общего средней зп за год') +
ylab('Логарифм среднего количества пропущенных дней за месяц')
ggplot(data, mapping = aes(x = monthly_income, y = sick_days, color = operations)) +
geom_point(alpha = 0.5) +
scale_x_log10() +
scale_y_log10() +
scale_color_discrete(name = "Наличие операций",labels=c("Да", "Нет")) +
ggtitle('Соотношение логарифма количества пропущенных рабочих дней\n на логарифм дохода',
subtitle = 'Распределение по наличию хирургических операций за последние 12 месяцев') +
xlab('Логарифм общего дохода за последний месяц') +
ylab('Логарифм среднего количества пропущенных дней за месяц')
ggplot(data = data)+
geom_bar(aes(x = gender, fill = gender))+
ggtitle('Распределение полов в выборке') +
scale_fill_discrete(name = "Пол",labels=c("Мужской", "Женский")) +
xlab('Пол') +
ylab('Количество')
ggplot(data = data)+
geom_bar(aes(x = status, fill = status))+
ggtitle('Распределение типов населенного пункта в выборке') +
scale_fill_discrete(name = "Тип",labels=c("Областной центр", "Город", "ПГТ", "Село")) +
xlab('Тип') +
ylab('Количество')
ggplot(data = data)+
geom_bar(aes(x = age))+
ggtitle('Распределение возрастов в выборке') +
xlab('Возраст') +
ylab('Количество')
ggplot(data = data)+
geom_bar(aes(x = children_count))+
ggtitle('Распределение по количеству детей младше 18-ти лет в выборке') +
xlab('Количество детей') +
ylab('Количество респондентов')
ggplot(data = data)+
geom_bar(aes(x = official_work, fill = official_work))+
ggtitle('Распределение количества работающих по трудовому договору в выборке') +
scale_fill_discrete(name = "Да/Нет",labels=c("Да", "Нет")) +
xlab('Наличие договора') +
ylab('Количество респондентов')
ggplot(data = data)+
geom_bar(aes(x = entrepreneurship, fill = entrepreneurship))+
ggtitle('Распределение количества занимающихся\n предпринимательством в выборке') +
scale_fill_discrete(name = "Да/Нет",labels=c("Да", "Нет")) +
xlab('Предприминательская ли деятельность') +
ylab('Количество респондентов')
ggplot(data = data)+
geom_bar(aes(x = disability, fill = disability))+
ggtitle('Распределение количества человек имеющих инвалидность в выборке') +
scale_fill_discrete(name = "Да/Нет",labels=c("Да", "Нет")) +
xlab('Наличие инвалидности') +
ylab('Количество респондентов')
ggplot(data = data)+
geom_bar(aes(x = branch, fill = branch))+
ggtitle('Распределение по отраслям в выборке')+
xlab('Номер отрасли') +
ylab('Количество респондентов') +
scale_fill_discrete(name = "Отрасль")#, labels=c('ЛЕГКАЯ, ПИЩЕВАЯ ПРОМЫШЛЕННОСТЬ', 'ГРАЖДАНСКОЕ МАШИНОСТРОЕНИЕ', 'ВОЕННО-ПРОМЫШЛЕННЫЙ КОМПЛЕКС', 'НЕФТЕГАЗОВАЯ ПРОМЫШЛЕННОСТЬ', 'ДРУГАЯ ОТРАСЛЬ ТЯЖЕЛОЙ ПРОМЫШЛЕННОСТИ', 'СТРОИТЕЛЬСТВО', 'ТРАНСПОРТ, СВЯЗЬ', 'СЕЛЬСКОЕ ХОЗЯЙСТВО', 'ОРГАНЫ УПРАВЛЕНИЯ', 'ОБРАЗОВАНИЕ', 'НАУКА, КУЛЬТУРА', 'ЗДРАВООХРАНЕНИЕ', 'АРМИЯ, МВД, ОРГАНЫ БЕЗОПАСНОСТИ', 'ТОРГОВЛЯ, БЫТОВОЕ ОБСЛУЖИВАНИЕ', 'ФИНАНСЫ', 'ЭНЕРГЕТИЧЕСКАЯ ПРОМЫШЛЕННОСТЬ', 'ЖИЛИЩНО-КОММУНАЛЬНОЕ ХОЗЯЙСТВО', 'ОПЕРАЦИИ С НЕДВИЖИМОСТЬЮ', 'СОЦИАЛЬНОЕ ОБСЛУЖИВАНИЕ', "ЮРИСПРУДЕНЦИЯ", 'ЦЕРКОВЬ', "ХИМИЧЕСКАЯ ПРОМЫШЛЕННОСТЬ", "ДЕРЕВООБРАБАТЫВАЮЩАЯ ПРОМЫШЛЕННОСТЬ, ЛЕСНОЕ ХОЗЯЙСТВО", "СПОРТ, ТУРИЗМ, РАЗВЛЕЧЕНИЯ", "УСЛУГИ НАСЕЛЕНИЮ", "IT, ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ", "ЭКОЛОГИЯ, ЗАЩИТА ОКРУЖАЮЩЕЙ СРЕДЫ", "ОРГАНИЗАЦИЯ ОБЩЕСТВЕННОГО ПИТАНИЯ", "СМИ, ИЗДАТЕЛЬСТВО, ПЕЧАТЬ, ТЕЛЕКОММУНИКАЦИИ", "РЕКЛАМА, МАРКЕТИНГ", "ОБЩЕСТВЕННЫЕ ОРГАНИЗАЦИИ, СОВЕТ ВЕТЕРАНОВ И ПР.")) +
model = lm(log(sick_days) ~ log(monthly_income), data)
summary(model)
##
## Call:
## lm(formula = log(sick_days) ~ log(monthly_income), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.44875 -0.47289 -0.05522 0.50974 3.06295
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.1674 0.3893 2.999 0.00276 **
## log(monthly_income) -0.1222 0.0376 -3.249 0.00119 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8335 on 1305 degrees of freedom
## Multiple R-squared: 0.008022, Adjusted R-squared: 0.007262
## F-statistic: 10.55 on 1 and 1305 DF, p-value: 0.001189