NIR 3

install.packages("devtools")
devtools::install_github("https://github.com/bdemeshev/rlms", force = TRUE)

library("memisc")
library("GGally")
library("dplyr")
library("psych")
library("lmtest")
library("sjPlot")
library("sgof")
library("ggplot2")
library("foreign")
library("car")
library("hexbin")
library("rlms")
library("devtools")
library("rstatix")
library("sandwich")
library("haven")



data <- read.csv("r17i_os26b.csv", sep=",", dec = ".", header=TRUE)
glimpse(data)
#Выборка данных для описания соц-эконом положения граждан РФ
#mh5 Пол(1-М, 2-Ж)
#mj13.2 Среднемесячная з/п за год
#m_marst Семейное положение
# (1-никогда не состоял в браке,
# 2-Состоит в зарег. браке,
# 3-Живут вместе, но не зарег,
# 4-Разведен и не состоит в браке,
# 5-Вдовец(вдова),
# 6-офиц. зарег, но живут не вместе)
#m_diplom Высшее образование
# (1-Окончил 0-6 классов,
# 2-Незаконч среднее образование(7-8кл),
# 3-2+что-то ещё,
# 4-Законч. Ср.Обр,
# 5-Законч. Ср.Спец.Обр,
# 6-Законч. Высш.Обр. и выше)
#m_age возраст
#status Тип населённого пункта
# (1-Обл.Центр,
# 2-Город,
# 3-Посёлок городского типа,
# 4-Село)
#mj6.2 Длительность рабочей недели

#отбираем нужные переменные
data = select(data, mh5, m_age, m_marst, m_diplom, status, mj13.2, mj6.2)
glimpse(data)
#убираем все  объекты с NA
data = na.omit(data)
#получаем представление наших даных(вывод)
glimpse(data)

#новая база данных с отфильтрованными значениями data2:
data2 = select(data)

#Cоздадим  дамми_перменные по параметру m_marst(семейное положение)
#1- женат, 0 - не женат
data2$wed1 = 0 #по умолчанию не женат
data2$wed1[which(data$m_marst == 2)] <- 1 #Если опрашиваемый женат присваеваем 1
data2$wed1[which(data$m_marst == 6)] <- 1 #Если опрашиваемый oфиц. зарег, но живут не вместе

#1 -  вдовец или разведен, 0 - не вдовец и не разведен
data2$wed2 = 0 #по умолчанию не вдовец и не разведен
data2$wed2[which(data$m_marst == 4)] <- 1 #Если респондент разведен и не состоит в браке присваеваем 1
data2$wed2[which(data$m_marst == 5)] <- 1 #Если респондент вдовец(вдова) присваеваем 1

#1 - никогда не состоял в браке, 0 - состоял
data2$wed3 = 0
data2$wed3[which(data$m_marst == 1)] <- 1 #Если респондент никогда не состоял в браке присваеваем 1

#проверим отсуствие линейной зависимости

vif(lm(data$m_marst ~ data2$wed1 + data2$wed2 + data2$wed3)) #все хорошо, зависимости нет

#делаем переменную пола(sex)
data2["sex"] = 0
data2$sex[which(data$mh5 == 1)] <- 1 #1 - мужчина, 0 - женщина

#делаеме из параметра status дамми-переменную, где 1- это город или областной центр и 0 в противоположном случае
data2$city_status = 0
data2$city_status[which(data$status == 1)] <- 1
data2$city_status[which(data$status == 2)] <- 1

#вводим параметр  higher_educ, характеризующий наличие полного высшего образования, где 1 - да, 0 - нет
data2$her_educ = 0
data2$her_educ[which(data$m_diplom == 2)] <- 1

#Нормализуем перемнные: зарплата, длительность рабочей недели и возраст

#зп
salary = data$mj13.2
data2$salary = (salary - mean(salary)) / sqrt(var(salary))

#длительность рабочей недели
wh = data$mj6.2
data2$wh = (wh - mean(wh)) / sqrt(var(wh))

#возраст
age = data$m_age
data2["age"] = (age - mean(age)) / sqrt(var(age))


glimpse(data2)
#(1)
model11 = lm(data = data2 ,salary ~ wh + age + city_status + her_educ + sex + wed1 + wed2 + wed3)
summary(model11)
vif(model11)
#Multiple R-squared:  0.03363
#Adjusted R-squared:  0.03178


#уберем  sex и city_status тк у них плохая  p-статистика(эта норм)
model12 = lm(data = data2 ,salary ~ wh + age  + her_educ  + wed1 + wed2 + wed3)

summary(model12)
vif(model12)
#Multiple R-squared:  0.03268
#Adjusted R-squared:  0.03129

#2
pow = 0.1
model21 = lm(data = data2, salary ~age  + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow))
summary(model21)
vif(model21)
#Multiple R-squared:  0.03268
#Adjusted R-squared:  0.03107

model22 = lm(data = data2, salary ~ wh + age  + her_educ  + wed1 + wed2 + wed3 + I(Mod(age)^pow))
summary(model22)
vif(model22)
#Multiple R-squared:  0.03372
#Adjusted R-squared:  0.0321


model23 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model23)
vif(model23)
#Multiple R-squared:  0.03372
#Adjusted R-squared:  0.0321

pow = 0.2
model24 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model24)
vif(model24)
#Multiple R-squared:  0.03385
#Adjusted R-squared:  0.03224

pow = 0.3
model25 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model25)
vif(model25)
#Multiple R-squared:  0.03397
#Adjusted R-squared:  0.03236

pow = 0.4
model26 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model26)
vif(model26)
#Multiple R-squared:  0.03408
#Adjusted R-squared:  0.03246

pow = 0.5
model27 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model27)
vif(model27)
#Multiple R-squared:  0.03
417
#Adjusted R-squared:  0.03256


pow = 0.6
model28 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model28)
vif(model28)
#Multiple R-squared:  0.03425
#Adjusted R-squared:  0.03264

pow = 0.7
model29 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model29)
vif(model29)
#Multiple R-squared:  0.03431
#Adjusted R-squared:  0.0327

pow = 0.8
model2_10 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model2_10)
vif(model2_10)
#Multiple R-squared:  0.03437
#Adjusted R-squared:  0.03275

pow = 0.9
model2_11 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model2_11)
vif(model2_11)
#Multiple R-squared:  0.0344
#Adjusted R-squared:  0.03279


pow = 1.0
model2_12 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model2_12)
vif(model2_12)
#Multiple R-squared:  0.03443
#Adjusted R-squared:  0.03282

pow = 1.1
model2_13 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model2_13)
vif(model2_13)
#Multiple R-squared:  0.03445
#Adjusted R-squared:  0.03284

pow = 1.2
model2_14 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model2_14)
vif(model2_14)
#Multiple R-squared:  0.03446
#Adjusted R-squared:  0.03284

pow = 1.3
model2_15 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model2_15)
vif(model2_15)
#Multiple R-squared:  0.03446
#Adjusted R-squared:  0.03285

pow = 1.4
model2_16 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model2_16)
vif(model2_16)
#Multiple R-squared:  0.03445
#Adjusted R-squared:  0.03284


pow = 1.5
model2_17 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model2_17)
vif(model2_17)
#Multiple R-squared:  0.03444
#Adjusted R-squared:  0.03283

pow = 1.6
model2_18 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model2_18)
vif(model2_18)
#Multiple R-squared:  0.03443
#Adjusted R-squared:  0.03281

pow = 1.7
model2_19 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model2_19)
vif(model2_19)
#Multiple R-squared:  0.0344
#Adjusted R-squared:  0.03279

pow = 1.8
model2_20 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model2_20)
vif(model2_20)
#Multiple R-squared:  0.03438
#Adjusted R-squared:  0.03277

pow = 1.9
model2_21 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model2_21)
vif(model2_21)
#Multiple R-squared:  0.03435
#Adjusted R-squared:  0.03274

pow = 2.0
model2_22 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model2_22)
vif(model2_22)
#Multiple R-squared:  0.03432
#Adjusted R-squared:  0.03271

#Произведения
model2_23 = m(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh) * Mod(age)))
summary(model2_23)
#Multiple R-squared:   0.03357
#Adjusted R-squared:    0.03196
vif(model2_23)
#попробуем убрать  wed1 тк vif(wed1) = 3
model2_24 = lm(data = data2, salary ~ wh + age  + her_educ + wed2 + wed3 + I(Mod(wh) * Mod(age)))
summary(model2_24)
#Multiple R-squared:   0.03357
#Adjusted R-squared:    0.03196
model2_25 = lm(data = data2, salary ~ wh + age +I(Mod(wh) * Mod(age)))
summary(model2_25)
#Multiple R-squared:   0.03097
#Adjusted R-squared:    0.03028
model2_26 = lm(data = data2, salary ~ wh + age + I(Mod(wh) * Mod(age)))
summary(model2_26)
#Multiple R-squared:   0.03097
#Adjusted R-squared:    0.03028

#Логарифмы
model2_27 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(log(Mod(wh))))
summary(model2_27)
vif(model2_27)
#Multiple R-squared:   0.03268
#Adjusted R-squared:    0.03129

model2_28 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(log(Mod(age))))
summary(model2_28)
vif(model2_28)
#Multiple R-squared:   0.01072
#Adjusted R-squared:    0.009305

model2_28 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(log(Mod(age) * Mod(wh))))
summary(model2_28)
vif(model2_28)
#Multiple R-squared:   0.0198
#Adjusted R-squared:    0.01839

model2_28 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(log(Mod(age))) + I(log(Mod(wh))))
summary(model2_28)
vif(model2_28)
#Multiple R-squared:   0.03358
#Adjusted R-squared:    0.03196


#3 Выберем лучшие модели
pow = 1.3
model2_29 = lm(data = data2, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model2_29)
vif(model2_29)
#Multiple R-squared:  0.03446
#Adjusted R-squared:  0.03285

#4 Молодые люди
#5
#С высшим образование, не из города - (her_educ = 1, city_status = 0)
#Разведенные женщины, без высшего образования - (wed2 = 1,her_edc = 0)

pow = 1.3
data3 = subset(data2,her_educ = 1 )
data3 = subset(data2,city_status = 0 )
model5 = lm(data = data3, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model5)
#Multiple R-squared:   0.03446
#Adjusted R-squared:    0.03285

pow = 1.3
data4 = subset(data2,wed2 = 1 )
data4 = subset(data2,her_educ = 0 )
model51 = lm(data = data4, salary ~ age + her_educ  + wed1 + wed2 + wed3 + I(Mod(wh)^pow) + I(Mod(age)^pow))
summary(model51)
#Multiple R-squared:   0.03446
#Adjusted R-squared:    0.03285