02-regression.Rmd

---
output: html_document
editor_options: 
  chunk_output_type: console
---
```{r purl = FALSE, cache = TRUE, include=FALSE}
knitr::opts_knit$set(global.par = TRUE)
knitr::opts_chunk$set(warning = FALSE, message = FALSE, collapse = TRUE, out.width = '100%')
library(tidyverse)
library(vroom)
```

# Лінійна регресія {#reg}
Модель лінійної регресії:
$$\hat{y} = \beta_0 + \sum \beta_j x_j$$
$\beta_0$ - вільний коефіцієнт (*bias, intercept*)

$\beta_j$ - ваговий коефіцієнт (*weights*)

$\beta_0,  \beta_1,  \dots, \beta_n$ - параметри

Інколи можуть скорочувати запис формули:
$$\hat{y} = \beta_0 + \left \langle \beta, x \right \rangle$$

$$\beta = \begin{pmatrix}
\beta_1 \\ 
\beta_2 \\ 
\dots \\ 
\beta_n

\end{pmatrix}$$

$$x = \begin{pmatrix}
x_1 \\ 
x_2 \\ 
\dots \\ 
x_n

\end{pmatrix}$$

## Коли такі моделі можна використовувати?

### Номінативні показники: one-hot encoding

Припустимо, що один з наших показників $x_j$ номінативний і він приймає значення з певної множини значень $x_j \in \left \{ c_1, c_2, \dots, c_m \right \}$

Основна мета one-hot encoding створити нові бінарні змінні з якими модель регресії зможе працювати:

$b_i(x) = \left [ f_i(x) = c_i \right ]$

Нотація Айверсона:
$$\left [ \models  \right ] = 1$$
$$\left [ ⊭  \right ] = 0$$

$b_1(x), b_2(x), \dots, b_m(x)$ - нові показники.

Як в такому випадку буде виглядати модель?

$$\hat{y} = \beta_0 + \beta_1\left [ b(x) = c_1 \right ] + \dots + \beta_m\left [ b(x) = c_m \right ] + \dots$$

### Бінаризація
В випадках, коли в нас існує нелінійна залежність між $y$ та $x$ є сенс використати бінеризацію.
```{r echo=FALSE}
tribble(~x, ~y,
        1, 2,
        2, 4,
        3, 8,
        4, 12,
        5, 11,
        6, 9,
        7, 8,
        8, 6,
        9, 5.5,
        10, 5) %>% 
  ggplot(aes(x, y)) +
  geom_smooth(se = FALSE)
```

Будуємо певну сітку значень ${t_1, t_2, \dots, t_m}$, тоді нові показники задамо як:

$$b_i(x) = [t_{i-1} < x_j \leqslant t_{i}], \;\; i = 1,\dots,m+1$$
```{r echo=FALSE}
tribble(~x, ~y,
        1, 2,
        2, 4,
        3, 8,
        4, 12,
        5, 11,
        6, 9,
        7, 8,
        8, 6,
        9, 5.5,
        10, 5) %>% 
  ggplot(aes(x, y)) +
  geom_smooth(se = FALSE) +
  geom_vline(xintercept = seq(1, 10, length.out = 12), color = "red")
```

Лінійна модель набуває вигляду:
$$\hat{y} = \beta_1\left [ t_{i-1} < x_j \leqslant t_{i} \right ] + \dots + \beta_m\left [ t_{m} < x_j \leqslant t_{m+1} \right ] + \dots$$
Межі інтервалів можна подавати, як перцентилі.

### Текстові дані: bag of words

## Похибки в задачах регресії
$$L(y, \hat{y})$$
### Квадратична функція похибок
$$L(y, \hat{y}) = (y - \hat{y})^2$$

$$MSE(\hat{y}, x) = \frac{1}{l}\sum(y - \hat{y}(x))^2$$

$$RMSE(\hat{y}, x) = \sqrt{\frac{1}{l}\sum(y - \hat{y}(x))^2}$$
$$R^2(\hat{y}, x) = 1 - \frac{\sum(y - \hat{y})^2}{\sum(y - \overline{y})^2}$$

### Абсолютна функція похибок
$$L(y, \hat{y}) = |y - \hat{y}|$$

$$MAE(\hat{y}, x) = \frac{1}{l}\sum|y - \hat{y}|$$
Позитивні сторони: стійка до викидів
```{r echo=FALSE}
tibble(
  y = c(1, 2, 3, 4, 5, 100, 7),
  fit_1 = c(2, 1, 2, 5, 6, 7, 6),
  fit_2 = c(4, 5, 6, 7, 8, 10, 10)
) %>% 
  mutate(abs_1 = abs(y - fit_1), 
         sqr_1 = (y - fit_1)^2,
         abs_2 = abs(y - fit_2), 
         sqr_2 = (y - fit_2)^2) %>% 
  knitr::kable()
```

```{r echo=FALSE}
tibble(
  y = c(1, 2, 3, 4, 5, 100, 7),
  fit_1 = c(2, 1, 2, 5, 6, 7, 6),
  fit_2 = c(4, 5, 6, 7, 8, 10, 10)
) %>% 
  mutate(abs_1 = abs(y - fit_1), 
         sqr_1 = (y - fit_1)^2,
         abs_2 = abs(y - fit_2), 
         sqr_2 = (y - fit_2)^2) %>% 
  select(4:7) %>% 
  map_dbl(mean) %>% 
  knitr::kable()
```

Негативні сторони: похідна не має інформації про близькість екстремуму + не має похідної в нулі.

### Huber loss
Поєднання квадратичної та абсолютної функції. Необхідно підбирати дельту. Не має другої похідної.
$$
L_{\delta}(y, \hat{y})\left\{\begin{matrix}
\frac{1}{2}(y - \hat{y})^2, & |y - \hat{y}| <  \delta \\ 
\delta (|y - \hat{y}| - \frac{1}{2}\delta) & |y - \hat{y}| \geqslant   \delta   
\end{matrix}\right.
$$

### Log-Cosh
Використовується гіперболічний косинус
$$
L_{\delta}(y, \hat{y}) = log (cosh(y - \hat{y}))
$$

### MSLE
Mean squared logarithmic error
$$
L(y, \hat{y}) = (log(\hat{y} + 1) - log(y + 1))^2
$$

### Відносні функції помилок

#### MAPE
$$
L_{\delta}(y, \hat{y}) = |\frac{y - \hat{y}}{y}|
$$

#### SMAPE
$$
L_{\delta}(y, \hat{y}) = \frac{|y - \hat{y}|}{(|y| + |\hat{y}|)/2}
$$

### Квантильна функції
Можна регулювати штраф за завищення і заниження похибок

Функцію помилок потрібно підбирати в залежності від задачі.

## Перенавчання
Нерідко в моделі машинного навчання стикаються з ситуацією *перенавчання* --- якість моделі на нових даних значно гірша ніж на навчальній вибірці. Тому важливо щоб наша модель вміла *узагалювати* свої результати на нові дані.

Для візуалізації цього ефекту проведемо симуляцію:
```{r sim}
library(tidyverse)
library(patchwork)

set.seed(1234)

df <- tibble(x = seq(1, 2, 0.05),
             y = cos(1.5 * pi * x) + rnorm(x, 0, 0.1))

plots <- map(c(1, 4, 15), function(d){
  ggplot(df, aes(x, y)) +
    geom_point() + 
    # geom_smooth(se = FALSE, color = "#2F6B57") + 
    geom_smooth(method = "lm", formula = y ~ poly(x, d), se = FALSE, color = "#2F6B57") +
    ggtitle(paste("Poly ", d))
})

plots[[1]] / plots[[2]] / plots[[3]]
```

Візуалізація демострує, що прості моделі мають недостатню точність, а складні моделі занадто добре підлаштовуються під вибірку, через що стають непридатними для подальшого використання.

Існує декілька варіантів виходу з ситуації перенавчання:

- регуляризація: штрафування моделі за складність

- крос-валідація: побудова низки моделей на підвибірках даних

- збільшення розмірності вибірки (про цей варіант часто забувають).

## Оцінювання якості моделей
Адекватна оцінка якості моделі грунтується на підході відкладеної вибірки: розмічені дані (дані, з відомими відповідями) розбиваються на дві частини: навчальну та тестову. На навчальній вибірці модель навчається а на тестовій перевіряється її якість. Якщо показник якості моделі на тестовій вибірці задовольняє наші потреби, можемо вважати, що модель знайшла певні закономірності в даних.

Але бувають випадки, коли якість моделі залежить від "формату" розбиття на підвибірки: ми можемо отримати різну якість моделі, якщо розіб'ємо дані за інших пропорцій або іншого початкового значення генератора випадкових величин. Вирішити таку проблему можна за допомогою крос-валідації. Дані розбиваються на $k$ блоків $X_1, X_2,\dots,X_k$ приблизно однакового розміру. Після чого будується $k$ моделей $\hat{y_1},\dots,\hat{y_k}$, при чому кожна $k$ модель навчається на всіх блоках окрім $k$. Після чого кожна модель оцінюється оцінюється по блоку який не приймав в навчанні, а результати усереднюються:

$$CV = \frac{1}{k}\sum{Q(\hat{y_i}, X_i)}$$

Як отримати фінальну модель для подальшого використання? Два варіанта:

- Навчаємо модель на всій вибірці даних, її параметри будуть підібрані на більшій кількості спостережень і можемо сподіватися, що якість моделі зросте.

- Будуємо композицію моделі з $\hat{y_1},\dots,\hat{y_k}$: наприклад, усереднення прогнозів всіх моделей, що може привести до підвищення стійкості моделі.
```{r}
library(tidymodels)

set.seed(1234)
n <- 1e3

ll_df <- tibble(
  x = runif(n, 0, 3),
  y = exp(1 + 0.75 * x + rnorm(n, sd = 0.5))
)

split <- initial_split(ll_df, prop = 0.8)
split

cars_train <- training(split)
cars_test <- testing(split)

cars_train
cars_test
```

```{r}
lm_spec <- linear_reg() %>%
  set_mode("regression") %>%
  set_engine("lm")

lm_fit <- lm_spec %>% 
  fit(y ~ x, data = cars_train)

augment(lm_fit, new_data = cars_train) %>%
  yardstick::rmse(truth = y, estimate = .pred)

augment(lm_fit, new_data = cars_test) %>%
  yardstick::rmse(truth = y, estimate = .pred)

predict(lm_fit, new_data = cars_test)

predict(lm_fit, new_data = cars_test, type = "conf_int")

bind_cols(
  predict(lm_fit, new_data = ll_df),
  ll_df
) %>%
  select(y, .pred)
```

```{r}
poly_tuned_rec <- recipe(y ~ x, data = cars_train) %>%
  step_poly(x, degree = tune())

poly_tuned_wf <- workflow() %>%
  add_recipe(poly_tuned_rec) %>%
  add_model(lm_spec)

cars_folds <- vfold_cv(cars_train, v = 10)

degree_grid <- grid_regular(degree(range = c(1, 10)), levels = 10)
degree_grid <- tibble(degree = seq(1, 10)) # same

tune_res <- tune_grid(
  object = poly_tuned_wf, 
  resamples = cars_folds, 
  grid = degree_grid
)

autoplot(tune_res)

collect_metrics(tune_res)

show_best(tune_res, metric = "rmse", n = 1)

best_degree <- show_best(tune_res, metric = "rmse", n = 1)

final_wf <- finalize_workflow(poly_tuned_wf, best_degree)

final_wf

final_fit <- fit(final_wf, cars_train)

final_fit

bind_cols(
  predict(final_fit, new_data = ll_df),
  ll_df
) %>%
  select(y, .pred)

augment(final_fit, new_data = cars_train) %>%
  yardstick::rmse(truth = y, estimate = .pred)

augment(final_fit, new_data = cars_test) %>%
  yardstick::rmse(truth = y, estimate = .pred)
```


## Градієнтний спуск
Пізніше...

## Приклад використання
```{r}
fm0 <- lm(y ~ x, ll_df)
X <- model.matrix(fm0)
f1 <- function(b) with(ll_df, sum(log(cosh(y - X %*% b))))
res <- optim(coef(fm0), f1, method = "BFGS")
res$par

f2 <- function(b) with(ll_df, mean(abs((y - X %*% b))))
res <- optim(coef(fm0), f2, method = "BFGS")
res$par
```


```{r}
f1 <- function(b) with(ll_df, sum(log(cosh(y - X %*% b))))
f2 <- function(b) with(ll_df, mean(abs((y - X %*% b))))
func <- list(f1, f2)
fm0 <- lm(y ~ x, ll_df)
X <- model.matrix(fm0)
param <- 0

for (i in 1:length(func)) {
  param[i] <- tibble(optim(coef(fm0), func[[i]], method = "BFGS")$par)
}

param <- param %>% 
  bind_rows(lm_fit$fit$coefficients) %>% 
  rename(Intercept = 1) %>% 
  mutate(loss = c("Log-Cosh", "MAE", "MSE"))

ll_df %>% 
  ggplot(aes(x, y)) +
  geom_point(alpha = .2) +
  geom_abline(data = param, aes(intercept = Intercept, slope = x, color = loss), size = 1)
```

bookdown::render_book("index.Rmd", output_dir = "docs")