Пример с одним признаком

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

Представим, что нам нужно предсказать заработную плату человека (\(Y\)). Хорошим признаком может быть сколько лет человек получал образование (\(X_1\)). В данном примере у нас есть всего один признак.

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

Давайте смоделируем данные для 50 людей. И посмотрим на график рассеяния, где по оси \(x\) отложим количество лет обучения, а по оси \(y\) доход.

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

ggplot(df, aes(x = education, y = income)) +
  geom_point()

Видно, что есть линейная зависимость. То есть нашу переменную \(y\) можно линейно выразить через переменную \(x\).

\[income = w_0 + w_1 education\] Напомним, как выглядит ошибка на одном наблюдении:

\[income_i - (w_0 + w_1 education_i)\]

На графике это размер пунктирной линии. Видно, что ошибка может быть положительной (недопрогнозировали) или отрицательной (перепрогнозировали).

df <- df %>%
        mutate(predict = 10000 + 5000 * education)

ggplot(df, aes(x = education, y = income)) +
  geom_point()+
  geom_abline(intercept = 10000, slope = 5000, color = 'blue') +
  geom_segment(aes(x = education, y = predict, 
                   xend = education, yend = income),
               linetype = 'dashed')

В качестве функционала возьмем \(MSE\). Он будет выглядеть как-то так:

\[MSE = \frac{1}{n} \sum_{i=1}^{n} (income_i - w_0 - w_1 education_i)^2\]

Модели отличаются друг от друга параметрами \(w_0\) и \(w_1\). Следовательно перед нами стоит задача минимизации функции потерь. Минимизировать эту функцию нужно по параметрам \(w_0\) и \(w_1\).

\[MSE = \frac{1}{n} \sum_{i=1}^{n} (income_i - w_0 - w_1 education_i)^2 \rightarrow \min_{w_0,w_1}\]

Возникает вопрос, какие значения дожны принимать \(w_0\) и \(w_1\)? На графике, это означает как именно нам нужно провести линию через точки, чтобы наша ошибка оказалась как можно меньше:

ggplot(df, aes(x = education, y = income)) +
  geom_point() + 
  geom_abline(intercept = 25000, slope = 5000, color = 'red') + 
  geom_abline(intercept = 10000, slope = 5000, color = 'blue') + 
  geom_abline(intercept = 9000, slope = 5500, color = 'green') + 
  geom_abline(intercept = 11000, slope = 4700, color = 'orange')

Обучим линейную регрессию.

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

В переменной model находится наша модель. Она имеет тип list. Из него можно получить значение коэффициентов или ошибки, совершенные моделью.

model$coefficients
## (Intercept)   education 
##   13964.261    4742.286
model$residuals
##           1           2           3           4           5           6 
##    155.9946  -6302.1452   6042.2882   1585.5488  -1229.5953   -490.3425 
##           7           8           9          10          11          12 
##   2425.5015  -1281.2759 -15538.3754   -152.8140  -2640.1377    696.5963 
##          13          14          15          16          17          18 
##    853.4759   5583.5768  -4908.2165   5705.8344    466.4363   2825.3994 
##          19          20          21          22          23          24 
##   2416.2114   3423.6414  -6166.4447  -2093.4952   2372.9261  -4729.2198 
##          25          26          27          28          29          30 
##  -2256.6154   2232.1723   3325.7287   1095.8660  -4437.5493  -6006.2176 
##          31          32          33          34          35          36 
##   7689.2182   1612.2046   1003.4002   -822.2898  -5546.9526   2407.7868 
##          37          38          39          40          41          42 
##  -1188.8614   -982.8221   3860.2974   4633.7736   7641.4219  -1923.9330 
##          43          44          45          46          47          48 
##   4353.9303   7360.7701  -5411.0399  -3164.3482  -5092.0070  -5498.4652 
##          49          50 
##   1299.1808   4793.9817

Также можно получить более подробную информацию с помощью функции summary.

summary(model)
## 
## Call:
## lm(formula = income ~ education, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15538.4  -2544.3    581.5   2725.4   7689.2 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  13964.3     2232.9   6.254 1.03e-07 ***
## education     4742.3      159.2  29.783  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4550 on 48 degrees of freedom
## Multiple R-squared:  0.9487, Adjusted R-squared:  0.9476 
## F-statistic:   887 on 1 and 48 DF,  p-value: < 2.2e-16

Интерпретация полученных коэффициентов

Из модели мы видим:

  • Cвободный коэффициент равен 13964.3. Это означает, что в среднем человек, который проучился 0 лет будет зарабатывать 13964.3 рубля.
  • Коэффициент перед переменной education равен 4742.3. Это означает, что если человек учится на один год больше, то в среднем он начинает зарабывать на 4742.3 рубля больше.

Эти оценки могут быть неточными. Чтобы получить доверительный интервалы, можно использовать функцию confint.

confint(model)
##                2.5 %    97.5 %
## (Intercept) 9474.648 18453.873
## education   4422.140  5062.432

Выбросы

Выбросы очень сильно влияют на линейную регрессию.

set.seed(42)
n <- 50
df <- tibble(education = seq(7, 20, length.out = n) + rnorm(n, sd = 2),
               income = 10000 + 5000 * education + rnorm(n, sd = 5000))
df <- add_row(df, education=c(19, 20, 21), income=c(1000, 1400, 2000))
ggplot(df, aes(x = education, y = income)) +
  geom_point()

На графике видно, что видимо кто-то указал зарплату не в рублях, а в долларах. Попробуем оценить линейную регрессию и нарисовать ее.

model <- lm(income ~ education, data = df)
model
## 
## Call:
## lm(formula = income ~ education, data = df)
## 
## Coefficients:
## (Intercept)    education  
##       37295         2611
coeff <- model$coefficients
ggplot(df, aes(x = education, y = income)) +
  geom_point() +
  geom_smooth(method = 'lm', se=F)

Видно, что линия регрессии проходит не через основное облако точек.

Использование большего количества признаков

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

n <- 100
df <- tibble(x1 = seq(-5, 5, length.out = n) + rnorm(n),
             x2 = seq(-10, 10, length.out = n) + rnorm(n),
             x3 = seq(12, 16, length.out = n) + rnorm(n),
             y = 5 + 2*x1 - 3*x2)

Для этого в формуле вам нужно просто перечислить те, которые вы хотите использовать.

model <- lm(y ~ x1 + x2 + x3, data = df)
summary(model)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3, data = df)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -4.255e-14 -1.313e-15 -3.270e-16  1.294e-15  3.416e-14 
## 
## Coefficients:
##               Estimate Std. Error    t value Pr(>|t|)    
## (Intercept)  5.000e+00  9.760e-15  5.123e+14   <2e-16 ***
## x1           2.000e+00  6.254e-16  3.198e+15   <2e-16 ***
## x2          -3.000e+00  3.402e-16 -8.819e+15   <2e-16 ***
## x3          -7.312e-17  6.935e-16 -1.050e-01    0.916    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.196e-15 on 96 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 1.286e+32 on 3 and 96 DF,  p-value: < 2.2e-16