Введение

В файле car_sales.csv содержится информация о продажах автомобилей различных марок:
- manufact - наименование производителя;
- model - название модели;
- country - страна производителя;
- sales - объем продаж в тысячах;
- resale - цена продажи 4-летнего автомобиля;
- type - тип (автомобиль или грузовик); - price - цена в тыс. долларов;
- engine_s - объем двигателя;
- horsepow - количество лошадиных сил;
- wheelbas - колесная база;
- width - ширина;
- length - длина;
- curb_wgt - вес кузова;
- fuel_cap - вместимость топливного бака в галлонах;
- mpg - эффективность расхода топлива, милль на галлон.

library(readr) # считывание данных из текстовых файлов
library(tidyverse) # манипулирование данными 
library(ggplot2) # визуализация 
library(ggfortify) # визуализация диагностических графиков
library(modelr) # вспомогательные функции для работы с моделями
library(broom) # преобразование результатов моделирования в табличный вид
library(GGally) # построение матрицы диаграмм рассеяния 
library(car) # функции для степенных преобразований
library(forcats) # работа с факторами

Sys.setlocale (category = "LC_ALL", locale = "pt_PT.UTF-8")
## [1] "pt_PT.UTF-8/pt_PT.UTF-8/pt_PT.UTF-8/C/pt_PT.UTF-8/C"

Задания

  1. Загрузите данные и обработайте пропуски, если они есть.
car_sales <- read.csv('data/car_sales.csv', skip = 20, sep = "\t", dec = ",")

car_sales %>%
  gather(key = 'Переменная', value = 'Значение') %>%
  group_by(`Переменная`) %>%
  summarise(`Пропуски` = sum(is.na(`Значение`))) %>%
  arrange(-`Пропуски`)
## Warning: attributes are not identical across measure variables;
## they will be dropped
## # A tibble: 15 x 2
##    Переменная Пропуски
##    <chr>         <int>
##  1 resale           36
##  2 mpg               3
##  3 curb_wgt          2
##  4 price             2
##  5 engine_s          1
##  6 fuel_cap          1
##  7 horsepow          1
##  8 length            1
##  9 wheelbas          1
## 10 width             1
## 11 country           0
## 12 manufact          0
## 13 model             0
## 14 sales             0
## 15 type              0
cars2 <- na.omit(car_sales)

Строчек, в который содержатся пропуски - 40. Это значительная цифра относительно общего числа наблюдений. Однако замена пропущенных значений средним, медианой или модой — грубый способ работы с ними. В данно ситуации, когда вариация данных невелика, такая грубая аппроксимация, возможно, приемлема, но если эта переменная влияет на выходную, это может исказить результаты.

В результате было принято решение удалить строки с рпопусками.

  1. С помощью матриц корреляции и визуализации, определите, какие факторы влияют на цену продажи 4-летнего автомобиля.

Отберем только числовые столбцы и исследуем корреляции между ними.

cars_num <- cars2 %>% select_if(is.numeric)
ggpairs(cars_num, 
        lower = list(continuous = wrap("smooth_lm", color = 'blue')))

# Коэффициенты корреляции и из значимость
cor(as.matrix(cars_num))
##                sales      resale       price    engine_s   horsepow
## sales     1.00000000 -0.27542560 -0.25170458  0.03811128 -0.1525377
## resale   -0.27542560  1.00000000  0.95475719  0.52718729  0.7731103
## price    -0.25170458  0.95475719  1.00000000  0.64917029  0.8534551
## engine_s  0.03811128  0.52718729  0.64917029  1.00000000  0.8616183
## horsepow -0.15253770  0.77311034  0.85345515  0.86161827  1.0000000
## wheelbas  0.40683932 -0.05368456  0.06704229  0.41001981  0.2259046
## width     0.17780194  0.17812845  0.30129205  0.67175604  0.5072748
## length    0.27233587  0.02539042  0.18259200  0.53734301  0.4009676
## curb_wgt  0.06718410  0.36327419  0.51139974  0.74283052  0.5986032
## fuel_cap  0.13804507  0.32479613  0.40649623  0.61686197  0.4797898
## mpg      -0.06671507 -0.39909720 -0.47992964 -0.72482938 -0.5959632
##             wheelbas      width      length   curb_wgt   fuel_cap
## sales     0.40683932  0.1778019  0.27233587  0.0671841  0.1380451
## resale   -0.05368456  0.1781284  0.02539042  0.3632742  0.3247961
## price     0.06704229  0.3012921  0.18259200  0.5113997  0.4064962
## engine_s  0.41001981  0.6717560  0.53734301  0.7428305  0.6168620
## horsepow  0.22590465  0.5072748  0.40096761  0.5986032  0.4797898
## wheelbas  1.00000000  0.6755589  0.85366861  0.6756090  0.6586544
## width     0.67555888  1.0000000  0.74322568  0.7359566  0.6721906
## length    0.85366861  0.7432257  1.00000000  0.6843051  0.5625042
## curb_wgt  0.67560895  0.7359566  0.68430507  1.0000000  0.8479935
## fuel_cap  0.65865437  0.6721906  0.56250422  0.8479935  1.0000000
## mpg      -0.47068533 -0.5996131 -0.46572235 -0.8189728 -0.8086330
##                  mpg
## sales    -0.06671507
## resale   -0.39909720
## price    -0.47992964
## engine_s -0.72482938
## horsepow -0.59596320
## wheelbas -0.47068533
## width    -0.59961305
## length   -0.46572235
## curb_wgt -0.81897279
## fuel_cap -0.80863299
## mpg       1.00000000

Наиболее значимыми факторами, влияющими на цену продажи 4-летнего автомобиля, представляются такие как resale, horsepower, engine_s.

  1. Выберите факторы и постройте несколько вариантов модели множественной регрессии. Влияет ли на характер изменения цены продажи 4-летнего транспортного средства тип транспортного средства? Страна производителя? Постройте модели с учетом этих факторов. Не забывайте проверять адекватность модели (значимость коэффициентов, выполнение допущений модели линейной регрессии). Запишите уравнения полученных моделей.

Линейная регрессия

На основе результатов корреляционного анализа можно предположить, что наилучшим предиктором для цены продажи 4-летнего автомобиля будет resale, т.к. коэффициент корреляции для этой переменной наибольший по модулю:

Построим модель простой линейной регрессии по переменной resale:

res_pr <- lm(resale ~ price, data = cars2)
tidy(res_pr) %>% select(term:std.error)
## # A tibble: 2 x 3
##   term        estimate std.error
##   <chr>          <dbl>     <dbl>
## 1 (Intercept)   -2.31     0.672 
## 2 price          0.783    0.0227
summary(res_pr)
## 
## Call:
## lm(formula = resale ~ price, data = cars2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.9015 -2.2156 -0.4313  2.0649 11.1464 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.30504    0.67202   -3.43 0.000839 ***
## price        0.78310    0.02275   34.43  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.466 on 115 degrees of freedom
## Multiple R-squared:  0.9116, Adjusted R-squared:  0.9108 
## F-statistic:  1185 on 1 and 115 DF,  p-value: < 2.2e-16

Получена модель: \(resale = -2.31+0.78*price\)

По F-критерию модель значима.

ggplot(add_predictions(cars2, res_pr, var = 'fit')) +
  geom_point(aes(price, resale)) +
  geom_line(aes(price, fit), colour = 'red') + 
  labs(title = 'Прогноз resale на основе price')

Подробная сводка по модели:

glance(res_pr) %>% select(r.squared, sigma)
## # A tibble: 1 x 2
##   r.squared sigma
## *     <dbl> <dbl>
## 1     0.912  3.47

Модель объяснила 91,2% дисперсии исходных данных, ее стандартная ошибка равна 3.47.

Множественная регрессия с двумя предикторами

Точность модели можно улучшить, если добавить в модель другие предикторы.

cars2 %>% 
  add_residuals(res_pr, var = "residual") %>% #добавили остатки
  select_if(is.numeric) %>% # убрали текстовые столбцы
  select(-resale, -price) %>% #убрали уже учтенные в модели переменные 
  gather('var', 'value', sales:fuel_cap) %>% #длинный формат данных
  group_by(var) %>% # сгруппировали по названию переменной
  summarise(r = cor(value, residual)) %>% # посчитали коэффициент корреляции для каждой группы
  arrange(desc(abs(r)))
## # A tibble: 8 x 2
##   var           r
##   <chr>     <dbl>
## 1 length   -0.501
## 2 curb_wgt -0.420
## 3 wheelbas -0.396
## 4 width    -0.368
## 5 engine_s -0.311
## 6 fuel_cap -0.213
## 7 horsepow -0.140
## 8 sales    -0.118

Наибольший по модулю коэффициент корреляции с остатками первой модели - у переменной length - длина машины.

Во второй модели рассмотрим одновременное влияние length и price.

pr_res_len <- lm(resale ~ price + length, data = cars2)
summary(pr_res_len)
## 
## Call:
## lm(formula = resale ~ price + length, data = cars2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.7409 -1.7684 -0.1832  2.0214  8.5962 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 21.33202    3.78468   5.636 1.27e-07 ***
## price        0.80617    0.01999  40.319  < 2e-16 ***
## length      -0.12911    0.02043  -6.320 5.23e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.996 on 114 degrees of freedom
## Multiple R-squared:  0.9345, Adjusted R-squared:  0.9334 
## F-statistic: 813.4 on 2 and 114 DF,  p-value: < 2.2e-16

Модель в целом значима. Получена модель:

\(resale = 21.33+0.81*price-0.13*length\)

Модель объяснила 93,4% дисперсии исходных данных вместо 91,2%.

ggplot(add_predictions(cars2, pr_res_len, var = 'fit'), 
       aes(fit, resale)) +
  geom_point() +
  geom_line(aes(fit, fit), color = 'red', linetype = 'dashed') +
  labs(title = 'Соответствие прогноза и факта для модели с length и price',
       x = 'Прогноз resale',
       y = 'Факт resale')

В реальности наблюдается отклонение фактических значений от предсказанных.

Регрессия с категориальными предикторами

Рассмотрим влияние 4-летней стоимости при перепродаже на цену. Поскольку в данных имеются два типа машин - автомобили и грузовики, можно предположить, что цена по-разному зависит для этих типов машин.

ggplot(cars2, aes(price, resale, colour = type)) +
  geom_point() +
  geom_smooth(method = lm, se = F) +
  labs(title = "Связь 4-летней стоимости при перепродаже и цены для двух типов машин",
       colour = "Тип машины")

Значения отличаются не кардинально, но для того, чтобы учесть в модели различие в типах машин, можно включить в модель фиктивную переменную, принимающую значение 0 или 1 в зависимости от типа машины. Значение 0 соответствует Automobile, а значение 1 - Truck. Умножение происходит из-за того, что угол наклона прямых различается.

cars2f <- cars2 %>%
  mutate(type = factor(type, levels = c('Automobile', 'Truck')))

res_type_factor <- lm(resale ~ price * type, data = cars2f)

summary(res_type_factor)
## 
## Call:
## lm(formula = resale ~ price * type, data = cars2f)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.319  -1.619  -0.228   2.176  10.727 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -2.64447    0.71077  -3.721 0.000311 ***
## price            0.80056    0.02308  34.694  < 2e-16 ***
## typeTruck        5.31165    2.13493   2.488 0.014306 *  
## price:typeTruck -0.23954    0.08270  -2.897 0.004531 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.367 on 113 degrees of freedom
## Multiple R-squared:  0.918,  Adjusted R-squared:  0.9158 
## F-statistic: 421.7 on 3 and 113 DF,  p-value: < 2.2e-16

Получена модель: \(resale = -2.64+0.80*price-0.24*price(truck=true)+5.31(truck=true)\)

Модель объяснила 91.8% дисперсии исходных данных.

Категориальные предикторы с более чем двумя категориями

Cтрана производителя машины имеет три уровня. В качестве базовой категории - Япония.

cars2f <- cars2f %>%
  mutate(country = factor(country, levels=c('Japan', 'Europe', 'USA')))

ggplot(cars2f, aes(price, resale, colour = country)) +
  geom_point() +
  geom_smooth(method = lm, se = F) + 
  labs(title = "Связь 4-летней стоимости при перепродаже и цены для машин из разных стран",
     colour = "Страна производителя")

Наклон прямых отличается слабо. Поэтому построим следующую модель:

res_country <- lm(resale ~ price + country, data = cars2f)
summary(res_country)
## 
## Call:
## lm(formula = resale ~ price + country, data = cars2f)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.3610 -1.5565 -0.2552  1.5150  9.8625 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.48584    0.74167  -0.655 0.513757    
## price          0.72808    0.02247  32.396  < 2e-16 ***
## countryEurope  3.70575    0.98707   3.754 0.000276 ***
## countryUSA    -1.67214    0.64064  -2.610 0.010277 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.055 on 113 degrees of freedom
## Multiple R-squared:  0.9325, Adjusted R-squared:  0.9307 
## F-statistic: 520.3 on 3 and 113 DF,  p-value: < 2.2e-16

Модель значима и объясняет 93.3% дисперсии исходных данных. Ее уравнение:

\(resale = -0.49+0.73*price+3.71(country = Europe)-1.67(country = USA)\)

  1. Сравните прогнозы, которые дают разные модели. Какая из моделей дает наиболее точные прогнозы?
glance(res_pr) %>% select(r.squared, sigma)
## # A tibble: 1 x 2
##   r.squared sigma
## *     <dbl> <dbl>
## 1     0.912  3.47
glance(pr_res_len) %>% select(r.squared, sigma)
## # A tibble: 1 x 2
##   r.squared sigma
## *     <dbl> <dbl>
## 1     0.935  3.00
glance(res_type_factor) %>% select(r.squared, sigma)
## # A tibble: 1 x 2
##   r.squared sigma
## *     <dbl> <dbl>
## 1     0.918  3.37
glance(res_country) %>% select(r.squared, sigma)
## # A tibble: 1 x 2
##   r.squared sigma
## *     <dbl> <dbl>
## 1     0.932  3.06

Таким образом, наиболее точный прогноз и наименьшая стандартная ошибка наблюдаются при множественной регрессии с двумя предикторами.

Исследуем остатки модели.

autoplot(pr_res_len)

Нелинейность выражена слабо. Наблюдаются выбросы на графиках. Наблюдения с наибольшими значениями остатков помечены номерами строк.Есть подозрение на гетероскедастичность и наличие влиятельных наблюдений.