Проект
Эконометрика
pacman::p_load(pacman, dplyr, ggplot2, tidyr, xtable, tidymodels,
psych, wooldridge, lmtest, sandwich, stargazer,
stats, modelsummary, readxl, GGally, ggthemes, viridis)auto.raw <- read.csv('auto_summary_ex.csv', stringsAsFactors = T)
auto <- auto.raw %>%
subset(auto.raw$cost > quantile(auto.raw$cost, 0.25) &
auto.raw$cost < quantile(auto.raw$cost, 0.75) + 1.5*IQR(auto.raw$cost))
auto$cost <- auto$cost/1000 # в тысячах рублей чтоб было
str(auto)## 'data.frame': 17015 obs. of 13 variables:
## $ X : int 5 6 10 18 19 24 31 40 44 48 ...
## $ mark : Factor w/ 95 levels "Acura","Alfa",..: 84 49 56 80 12 35 76 48 84 9 ...
## $ model : Factor w/ 3002 levels "(DongFeng-Honda) M-NV",..: 2120 2322 1156 1302 811 1046 2218 2274 2613 1221 ...
## $ mileage : int 98000 87000 88000 185000 128000 68000 90000 114000 62000 118000 ...
## $ year : int 2018 2017 2019 2011 2014 2013 2017 2017 2017 2016 ...
## $ litres : num 1.6 2 2 2 3 2.4 1.6 2 2 6.2 ...
## $ power : int 110 238 211 147 249 190 110 180 150 409 ...
## $ engine : Factor w/ 5 levels "Бензин","Газ",..: 1 1 3 4 1 1 1 4 4 1 ...
## $ transmission: Factor w/ 4 levels "автомат","вариатор",..: 3 1 1 3 1 1 3 1 4 1 ...
## $ body : Factor w/ 17 levels "внедорожник",..: 12 1 12 1 1 1 7 1 1 1 ...
## $ drive : Factor w/ 3 levels "задний","передний",..: 2 3 1 3 3 3 2 3 3 3 ...
## $ colour : Factor w/ 65 levels "Active (Пакет 1)",..: 57 57 63 56 63 38 58 56 56 63 ...
## $ cost : num 1230 3060 3740 1340 1230 2420 1140 2840 3010 4710 ...auto.raw %>% ggplot(aes(cost)) +
geom_boxplot(fill = 'cadetblue3') +
xlab('Стоимость автомобиля')
#ggsave('boxplot.png', height = 6, width = 20)У переменной litres максимум 795, а медиана 2. Убираем выбросы. Очевидно, объём двигателя чаще всего колеблется в пределах 2-5 литров. Но учитывая вариацию автомобилей оставим всё, что ниже 7.
print('До удаления выбросов')## [1] "До удаления выбросов"auto.raw$litres %>% describe## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 25498 3.13 23.28 1.8 1.9 0.3 0.7 1020 1019.3 26.86 806.97 0.15auto <- auto %>% subset(litres < 7)
print('После удаления выбросов')## [1] "После удаления выбросов"auto$litres %>% describe## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 16973 2.05 0.78 2 1.92 0.59 0.7 6.7 6 2.31 6.74 0.01auto %>%
ggplot(aes(litres))+
geom_boxplot()
datasummary(('Мощность двигателя' = power) + (Стоимость = cost) + ('Объём двигателя' = litres) + (Пробег = mileage) ~ Mean + SD + Min + Max + Median,
data = auto,
output = 'html',
#output = 'latex'
)| Mean | SD | Min | Max | Median | |
|---|---|---|---|---|---|
| “Мощность двигателя” | 160.34 | 72.22 | 29 | 760 | 144 |
| Стоимость | 1849.97 | 991.55 | 720.00 | 4830.00 | 1540.00 |
| “Объём двигателя” | 2.05 | 0.78 | 0.70 | 6.70 | 2.00 |
| Пробег | 104521.59 | 75032.03 | 0 | 692000 | 94000 |
datasummary(('Мощность двигателя' = power) + (Стоимость = cost) + ('Объём двигателя' = litres) + (Пробег = mileage) ~ Var + Mean + SD + Min + Max + Median,
data = auto.raw,
output = 'html'
#output = 'latex'
)| Var | Mean | SD | Min | Max | Median | |
|---|---|---|---|---|---|---|
| “Мощность двигателя” | 7115.66 | 158.01 | 84.35 | 13 | 900 | 136.00 |
| Стоимость | 8e+12 | 2064155.84 | 2850503.58 | 30000 | 55590000 | 1270000.00 |
| “Объём двигателя” | 541.74 | 3.13 | 23.28 | 0.70 | 1020.00 | 1.80 |
| Пробег | 7e+09 | 116547.49 | 84077.02 | 0 | 1170000 | 106000.00 |
Делаем дамми переменные для трансмиссии
dummies.auto <- recipe(~., data = auto) %>%
step_dummy(c(transmission), one_hot = T) %>%
prep(training = auto) %>%
bake(new_data = NULL)
dummies.auto$OtherTr <- dummies.auto$transmission_вариатор + dummies.auto$transmission_робот
dummies.auto <- dummies.auto[ , -c(14, 16)]
dummies.auto <- dummies.auto %>% rename('Automatic' = transmission_автомат,
'Manual' = transmission_механика)
dummies.auto[, 13:15]## # A tibble: 16,973 × 3
## Automatic Manual OtherTr
## <dbl> <dbl> <dbl>
## 1 0 1 0
## 2 1 0 0
## 3 1 0 0
## 4 0 1 0
## 5 1 0 0
## 6 1 0 0
## 7 0 1 0
## 8 1 0 0
## 9 0 0 1
## 10 1 0 0
## # … with 16,963 more rowsСколько разных классов?
table(auto$transmission)##
## автомат вариатор механика робот
## 7191 1475 6545 1762
linreg_0 <- lm(cost ~ power, data = dummies.auto)
cov1 <- vcovHC(linreg_0, type = 'HC0')
nelinreg <- lm(cost ~ power + power:Automatic + power:Manual,
data = dummies.auto)
cov2 <- vcovHC(nelinreg, type = 'HC0')
nelinreg.contr <- lm(cost ~ power + drive + engine + litres + mileage,
data = dummies.auto)
cov3 <- vcovHC(nelinreg.contr, type = 'HC0')
super.nelinreg.contr <- lm(cost ~ power + power:Automatic + power:Manual + drive + engine + litres + mileage,
data = dummies.auto)
cov4 <- vcovHC(super.nelinreg.contr, type = 'HC0')Предельные эффекты
msr <- margins::margins_summary(super.nelinreg.contr,
variables = 'power',
at = list(Automatic = c(0, 1),
Manual = c(0, 1),
OtherTr = c(0, 1)), # сравнение по трансмиссии
vcov = vcovHC(super.nelinreg.contr, type = 'HC0'))[c(3,5),]
msr%>%
ggplot() +
geom_point(aes(Manual, AME), size = 3) +
geom_errorbar(aes(Manual, AME, ymin = lower, ymax = upper), width = .1, linewidth = 2) +
scale_x_discrete(
breaks = c('0', '1'),
labels = c('Automatic', 'Manual')
) +
xlab('Transmission') +
geom_text(x = 1.2, y = 3.13,
label = "Manual", size = 9) +
geom_text(x = 0.25, y = 4.3,
label = "Automatic", size = 9) +
theme(legend.position="none")
#ggsave('AME.png', width = 12, height = 7)msr[, c(2, 3, 5, 8, 9, 10)]## Automatic Manual AME
## 0 1 3.109
## 1 0 4.329
## p
## 0.00000000000000000000000000000019647690284662775204310714016742167586926370859
## 0.00000000000000000000000000000000000000000000000000000000000000000000000006769
## lower upper
## 2.586 3.631
## 3.863 4.796modelsummary(models = list(linreg_0,
nelinreg.contr,
nelinreg,
super.nelinreg.contr),
vcov = list(cov1,
cov3,
cov2,
cov4),
#output = 'latex',
#output = 'text',
coef_omit = "mileage|driveполный|driveпередний|engineГаз|engineГибрид|engineДизель|engineЭлектро|litres|power:transmissionвариатор|power:transmissionмеханика|power:transmissionробот",
coef_rename = c("power" = "Можность двигателя в лошадиных силах",
"driveполный" = "Полный привод",
"driveпередний" = "Передний привод",
"engineГаз" = "Газовый двигатель",
"engineГибрид" = "Гибридный двигатель",
"engineДизель" = "Дизельный двигатель",
"engineЭлектро" = "Электро двигатель",
"litres" = "Объём двигателя",
"power:transmissionвариатор" = "Мощность*Вариатор",
"power:transmissionмеханика" = "Мощность*Механика",
"power:transmissionробот" = "Мощность*Робот",
'(Intercept)' = 'Константа',
'power:Automatic' = 'Мощность*Автомат',
'power:Manual' = 'Мощность*Механика',
'mileage' = 'Пробег'),
statistic = "std.error",
stars = TRUE,
gof_omit = ".*",
notes = list("В скобках даны робастные стандартные ошибки",
"Все регрессии содержат контрольные переменные: Mileage (пробег), Litres (объём двигателя)\n, Engine (тип двигателя), Drive (тип привода)"),
title = "Результаты оценивания")| Model 1 | Model 2 | Model 3 | Model 4 | |
|---|---|---|---|---|
| Константа | 1068.524*** | 1582.174*** | 1188.241*** | 1672.305*** |
| (19.342) | (44.315) | (22.227) | (45.176) | |
| Можность двигателя в лошадиных силах | 4.874*** | 4.615*** | 4.910*** | 4.399*** |
| (0.122) | (0.236) | (0.177) | (0.250) | |
| Мощность*Автомат | -0.410** | -0.070 | ||
| (0.129) | (0.123) | |||
| Мощность*Механика | -1.917*** | -1.290*** | ||
| (0.136) | (0.130) | |||
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||||
| В скобках даны робастные стандартные ошибки | ||||
| Все регрессии содержат контрольные переменные: Mileage (пробег), Litres (объём двигателя) , Engine (тип двигателя), Drive (тип привода) |
pow.std <- sqrt(cov4[2, 2]) # SE(beta_1_hat)
pm.std <- sqrt(cov4[11, 11]) # SE(beta_2_hat)
pa.std <- sqrt(cov4[10, 10]) # SE(beta_3_hat)
cov.am <- cov4[10, 11]
beta_1 <- super.nelinreg.contr$coefficients[2]
beta_2 <- super.nelinreg.contr$coefficients[10]
beta_3 <- super.nelinreg.contr$coefficients[11]Первая гипотеза \(\beta_1 > 0\). Применим односторонний z-тест. Для этого посчитаем t - статистику = \(\dfrac{\hat\beta_1 - \beta_1}{SE(\hat\beta_1)}\)
t.stat.pow <- beta_1/pow.std
pval1 <- pnorm(q=t.stat.pow,
lower.tail = F)
cat('t-статистика = ', t.stat.pow, '\np-value = ', round(pval1))## t-статистика = 17.57431
## p-value = 0Вторая гипотеза
lh <- car::linearHypothesis(super.nelinreg.contr,
c("power = 0",
"power:Automatic = 0",
"power:Manual = 0"),
df = Inf,
test = "Chisq",
white.adjust = "hc0")
cat('значение статистики = ', lh$Chisq[2], '\np-value = ', round(lh$`Pr(>Chisq)`[2]))## значение статистики = 549.5699
## p-value = 0Третья гипотеза
std.sum <- pm.std + pa.std - 2*cov.am
t.stat.am <- (beta_2 - beta_3)/std.sum
pval2 <- pnorm(t.stat.am,
lower.tail = F)
cat('t-статистика = ', t.stat.am, '\np-value = ', round(pval2))## t-статистика = 5.24171
## p-value = 0