Линейная регрессия: часть 3

---

library(tibble)
library(dplyr)
library(tidyr)
library(ggplot2)

Прогноз

Когда мы построили модель, мы можем использовать ее для предсказания. Это можно сделать с помощью функции predict.

set.seed(42)
n <- 50
data <- tibble(education = seq(7, 20, length.out = n) + rnorm(n, sd = 2),
               income = 10000 + 5000 * education + rnorm(n, sd = 5000))

model <- lm(income ~ education, data = data)
model

## 
## Call:
## lm(formula = income ~ education, data = data)
## 
## Coefficients:
## (Intercept)    education  
##       13964         4742

data$pred_income <- predict(model, data)
data$pred_income

##  [1]  60163.22  43062.50  53120.69  56937.17  56027.20  52444.50  69045.35
##  [8]  55069.56  76369.41  57888.86  72117.97  82687.85  49085.40  60872.12
## [15]  63509.97  72064.34  64594.76  43353.60  46660.31  83585.96  69415.08
## [22]  56686.62  73209.16  87618.55  95331.14  74531.38  77432.27  64407.67
## [29]  86752.50  77576.75  89224.74  92848.22  97238.82  82904.05  94726.90
## [36]  74910.68  85013.67  85641.60  72072.52  96571.01  99440.37  95320.25
## [43] 107193.73  94368.55  89541.63 107882.44  97339.79 119990.35 103459.74
## [50] 115028.52

Графическая иллюстрация модели

Когда признаков много, невозможно представить себе многомерное пространство. Нарисовать это тоже невозможно. Но можно посмотреть на то, как работает модель. Попробуем отложить по оси \(X\) наш прогноз, а по оси \(Y\) реальное значение. Идеально спрогнозированные наблюдения лежат на прямой \(Y = X\).

ggplot(data, aes(x = pred_income, y = income)) +
  geom_point() +
  geom_abline(color = "blue") + 
  xlab('Прогноз') + 
  ylab('Истинное значение')

На этом графике можно видеть недопрогноз и перепрогноз.

Метрики качества

Чтобы оценить насколько хорошо отработала модель нужна какая-то мера качества модели. Для задачи регрессии основными являются:

1. RMSE

\[RMSE = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - \hat y_i)^2}\]

2. MAE

\[MAE = \frac{1}{n} \sum_{i=1}^{n} |y_i - \hat y_i|\]

3. MAPE

\[MAPE = \frac{1}{n} \sum_{i=1}^{n} |\frac{y_i - \hat y_i}{y_i}|\]

library(Metrics)

rmse(data$income, data$pred_income)

## [1] 4458.386

mae(data$income, data$pred_income)

## [1] 3514.527

Чтобы понять много это или мало нужно понять чему равно среднее значение таргета.

mean(data$income)

## [1] 77646.79

mape(data$income, data$pred_income)

## [1] 0.04743514

Нелинейная зависимость

Не всегда зависимости линейны.

n <- 100
df <- tibble(x = seq(-1, 5, length.out = n) + rnorm(n, sd = 0.1),
             y = x^2 - 2*x + rnorm(n))

ggplot(df, aes(x, y)) + 
  geom_point()

Строить линейную регрессию на таких данных не кажется разумным.

ggplot(df, aes(x, y)) + 
  geom_point() +
  geom_smooth(method = 'lm', se = F)

Красная линия более удачная модель для наших данных.

ggplot(df, aes(x, y)) + 
  geom_point() +
  geom_smooth(method = 'lm', se = F) + 
  geom_smooth(se = F, color = 'red')

Чтобы построить такую модель можно использовать полиномы признака x: квадрата, кубы и так далее. Давайте попробуем взять квадрат этой модели и обучить нашу модель.

Вариант 1

df <- df %>%
  mutate(x2 = x**2)
model <- lm(y ~ x + x2, data = df)
model

## 
## Call:
## lm(formula = y ~ x + x2, data = df)
## 
## Coefficients:
## (Intercept)            x           x2  
##    -0.09175     -1.86596      0.97371

Вариант 2

model <- lm(y ~ x + I(x^2), data = df)
model

## 
## Call:
## lm(formula = y ~ x + I(x^2), data = df)
## 
## Coefficients:
## (Intercept)            x       I(x^2)  
##    -0.09175     -1.86596      0.97371

Вариант 3

model <- lm(y ~ poly(x, 2, raw=T), data = df)
model

## 
## Call:
## lm(formula = y ~ poly(x, 2, raw = T), data = df)
## 
## Coefficients:
##          (Intercept)  poly(x, 2, raw = T)1  poly(x, 2, raw = T)2  
##             -0.09175              -1.86596               0.97371

Визуализируем нашу модель.

ggplot(df, aes(x, y)) + 
  geom_point() +
  geom_smooth(method = 'lm', se = F) +
  geom_smooth(method = 'lm', formula = y ~ poly(x, 2, raw=T), se = F, col='red')

Изначально мы не знали сколько полиномов нужно взять. Давайте возьмем больше.

ggplot(df, aes(x, y)) + 
  geom_point() +
  geom_smooth(method = 'lm', se = F) +
  geom_smooth(method = 'lm', formula = y ~ poly(x, 10, raw=T), se = F, col='red')

model <- lm(y ~ poly(x, 10, raw=T), data = df)

Видно, что такая модель начинает подгоняться под наши данные. В машинном обучении этот факт называется переобученим.

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

Train and test

Чтобы контролировать эффект переобучения принято делить данные на тренировочные(train) и тестовые(test). На тренировочном датасете мы учимся, а на тестовом смотрим результаты нашей модели и оцениваем насколько модель хороша. Обычно набор данных делят 70 (train) на 30 (test).

Сгененрируем очень простой датасет.

n <- 100
df <- tibble(x = seq(-3, 1, length.out = n) + rnorm(n, sd = 0.1),
             y = x^2 + rnorm(n, sd = 0.3))

ggplot(df, aes(x, y)) + 
  geom_point()

Попробуем поделить наши данные на 2 части.

df <- df %>%
  arrange(x) %>%
  mutate(type = c(rep('Train', n*0.7), rep('Test', n*0.3)))

Визуализируем.

ggplot(df, aes(x, y, color = type)) +
  geom_point()

Если обучить линейную регрессию на тренировочных данных, то мы получим следующий результат.

train <- df %>% filter(type == 'Train')
test <- df %>% filter(type == 'Test')

model <- lm(y ~ x, data = train)
coeff <- model$coefficients

ggplot(df, aes(x, y, color = type)) +
  geom_point() + 
  geom_abline(intercept = coeff[1], slope = coeff[2], color = 'darkgreen') + 
  ylim(-4, 10)

Прогноз на тестовых данных будет очень плохим. Это происходит, так как модель не видела данные при х больше 0. При делении данных на train и test нужно перемешивать ваши данные.

set.seed(42)
n <- 200
df <- tibble(x = seq(-3, 1, length.out = n) + rnorm(n, sd = 0.1),
             y = x^2 + rnorm(n, sd = 0.3))


library(rsample)

spliter <- initial_split(df, prop = 0.7)

train <- training(spliter)
test <- testing(spliter)

train$type <- 'Train'
test$type <- 'Test'

data <- bind_rows(train, test)

Визуализируем.

ggplot(data, aes(x, y, color = type)) +
  geom_point() +
  geom_smooth(data = train, se = F, formula = y ~ poly(x, 2, raw = T))