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

---

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

Функция lm

На прошлой посиделке мы узнали, что можно обучить линейную регрессию с помощью функции lm. Наши данные находятся в датасете data. Мы хотим предсказать заработную плату человека (income). Нашим единственныи признаком является количество лет, которые человек отучился (education).

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))

head(df)

##                                               
## 1 function (x, df1, df2, ncp, log = FALSE)    
## 2 {                                           
## 3     if (missing(ncp))                       
## 4         .Call(C_df, x, df1, df2, log)       
## 5     else .Call(C_dnf, x, df1, df2, ncp, log)
## 6 }

Нарисуем график рассеяния этих двух переменных. На этом графике мы видим линейную зависимость переменных.

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

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

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

В переменной 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 = data)
## 
## 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вободный коэффициент равен 8683. Это означает, что в среднем человек, который проучился 0 лет будет зарабатывать 8683 рубля.
• Коэффициент перед переменной education равен 5091. Это означает, что если человек учится на один год больше, то в среднем он начинает зарабывать на 5091 рубль больше.

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

confint(model)

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

Выбросы

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

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))
data <- add_row(data, education=c(19, 20, 21), income=c(1000, 1400, 2000))

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

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

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

## 
## Call:
## lm(formula = income ~ education, data = data)
## 
## Coefficients:
## (Intercept)    education  
##       37295         2611

coeff <- model$coefficients

ggplot(data, 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

Категориальные переменные

Давайте рассмотрим следующий набор данных.

n <- 50
data <- tibble(education = seq(7, 20, length.out = n) + rnorm(n, sd = 2),
               type = rep(c('college', 'university'), times = n/2),
               income = ifelse(type=='college', 
                               10000 + 5000 * education + rnorm(n, sd = 5000),
                               40000 + 5000 * education + rnorm(n, sd = 5000)))

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

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

Но если попробовать покрасить точки, используя переменную type, то можно понять, что доход зависит от типа образования, которое получил человек.

ggplot(data, aes(x = education, y = income, color = type)) +
  geom_point()

Бинарная переменная

Чтобы отобразить это в моделе можно закодировать нашу переменную type в 0, если образование college, и в 1, если образование university. Это можно сделать с помощью функции ifelse.

data <- data %>%    
  mutate(type_encode = ifelse(type == "college", 0, 1))

head(data)

## # A tibble: 6 x 4
##   education type       income type_encode
##       <dbl> <chr>       <dbl>       <dbl>
## 1      9.67 college    51830.           0
## 2      5.53 university 72208.           1
## 3      7.64 college    49063.           0
## 4      7.89 university 80150.           1
## 5      6.90 college    45046.           0
## 6      6.33 university 70655.           1

Попробуем обучить эту модель используя бинарную переменную.

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

## 
## Call:
## lm(formula = income ~ education + type_encode, data = data)
## 
## Coefficients:
## (Intercept)    education  type_encode  
##       12024         4885        29895

coeff <- model$coefficients

Коэффициент перед type_encode показывает на сколько в среднем человек с высшим образованием зарабатывает больше чем человек со средним образованием.

ggplot(data, aes(x = education, y = income, color = type)) +
  geom_point() +
  geom_abline(intercept = coeff[1], slope = coeff[2], color='darkred') +
  geom_abline(intercept = coeff[1] + coeff[3], slope = coeff[2], color='darkblue')

Label Encoding

Что делать, если переменная принимает больше двух значений?

n <- 60
data <- tibble(education = seq(7, 20, length.out = n) + rnorm(n, sd = 2),
               type = rep(c('college', 'university', 'PhD'), times = n/3),
               income = ifelse(type=='college', 
                               10000 + 5000 * education + rnorm(n, sd = 5000),
                               ifelse(type=='university',
                                      40000 + 5000 * education + rnorm(n, sd = 5000),
                                      60000 + 5000 * education + rnorm(n, sd = 5000))))

head(data)

## # A tibble: 6 x 3
##   education type        income
##       <dbl> <chr>        <dbl>
## 1      3.88 college     20058.
## 2      7.87 university  81797.
## 3      7.13 PhD         97475.
## 4      9.42 college     58450.
## 5      9.38 university  90140.
## 6      8.70 PhD        107226.

Визуализируем данные.

ggplot(data, aes(x = education, y = income, color = type)) +
  geom_point() 

Давайте просто закодируем переменную как в прошлый раз.

data <-  data %>%
    mutate(type_encode = 
            case_when(type == 'college' ~ 1,
                      type == 'university' ~ 2,
                      type == 'PhD' ~ 3))

head(data)

## # A tibble: 6 x 4
##   education type        income type_encode
##       <dbl> <chr>        <dbl>       <dbl>
## 1      3.88 college     20058.           1
## 2      7.87 university  81797.           2
## 3      7.13 PhD         97475.           3
## 4      9.42 college     58450.           1
## 5      9.38 university  90140.           2
## 6      8.70 PhD        107226.           3

Обучим модель.

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

## 
## Call:
## lm(formula = income ~ education + type_encode, data = data)
## 
## Coefficients:
## (Intercept)    education  type_encode  
##      -13907         5015        25022

coeff <- model$coefficients

Смотрим на результат.

ggplot(data, aes(x = education, y = income, color = type)) +
  geom_point() +
  geom_abline(intercept = coeff[1] + coeff[3]*1, slope = coeff[2], color='darkred') +
  geom_abline(intercept = coeff[1] + coeff[3]*2, slope = coeff[2], color='darkblue') +
  geom_abline(intercept = coeff[1] + coeff[3]*3, slope = coeff[2], color='darkgreen')

Такой подход не очень хорош. Попробуем закодировать наши переменные по-другому.

data <-  data %>%
    mutate(type_encode = 
            case_when(type == 'college' ~ 1,
                      type == 'university' ~ 3,
                      type == 'PhD' ~ 2))

head(data)

## # A tibble: 6 x 4
##   education type        income type_encode
##       <dbl> <chr>        <dbl>       <dbl>
## 1      3.88 college     20058.           1
## 2      7.87 university  81797.           3
## 3      7.13 PhD         97475.           2
## 4      9.42 college     58450.           1
## 5      9.38 university  90140.           3
## 6      8.70 PhD        107226.           2

Обучим модель.

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

## 
## Call:
## lm(formula = income ~ education + type_encode, data = data)
## 
## Coefficients:
## (Intercept)    education  type_encode  
##        7221         4945        14922

coeff <- model$coefficients

Смотрим на результат.

ggplot(data, aes(x = education, y = income, color = type)) +
  geom_point() +
  geom_abline(intercept = coeff[1] + coeff[3]*1, slope = coeff[2], color='darkred') +
  geom_abline(intercept = coeff[1] + coeff[3]*2, slope = coeff[2], color='darkblue') +
  geom_abline(intercept = coeff[1] + coeff[3]*3, slope = coeff[2], color='darkgreen')

В данном случае кодировка очень важна, потому что 1 и 2 ближе чем 1 и 3.

One Hot Encoding

Чтобы избежать проблем с Label Encoding можно использовать One Hot Encoding (OHE).

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

n <- 60
data <- tibble(education = seq(7, 20, length.out = n) + rnorm(n, sd = 2),
               type = rep(c('college', 'university', 'PhD'), times = n/3),
               income = ifelse(type=='college', 
                               10000 + 5000 * education + rnorm(n, sd = 5000),
                               ifelse(type=='university',
                                      40000 + 5000 * education + rnorm(n, sd = 5000),
                                      60000 + 5000 * education + rnorm(n, sd = 5000))))

head(data)

## # A tibble: 6 x 3
##   education type        income
##       <dbl> <chr>        <dbl>
## 1      7.52 college     42433.
## 2      7.54 university  78656.
## 3      9.30 PhD        109513.
## 4      7.54 college     42139.
## 5      7.98 university  85471.
## 6      5.96 PhD         92444.

Попробуем на нашем датасете.

library(onehot)

OHE <- function(df){
  encoder <- onehot(df)
  df_OHE <- as_tibble(predict(encoder, df))
  return(df_OHE)
}

data$type <- factor(data$type)
data_OHE <- OHE(data)

head(data_OHE)

## # A tibble: 6 x 5
##   education `type=college` `type=PhD` `type=university`  income
##       <dbl>          <dbl>      <dbl>             <dbl>   <dbl>
## 1      7.52              1          0                 0  42433.
## 2      7.54              0          0                 1  78656.
## 3      9.30              0          1                 0 109513.
## 4      7.54              1          0                 0  42139.
## 5      7.98              0          0                 1  85471.
## 6      5.96              0          1                 0  92444.

Переменных уже много и писать их вручную не очень. Вместо этого можно поставить точку. Тогда все переменные, которые есть в вашем наборе данных будут использованы. В данном случае мы хотим оценить нашу модель без свободного коэффициента, поэтому поставим -1.

model <- lm(income ~ . - 1 , data = data)
model

## 
## Call:
## lm(formula = income ~ . - 1, data = data)
## 
## Coefficients:
##      education     typecollege         typePhD  typeuniversity  
##           5091            9329           59416           38955

coeff <- model$coefficients

Смотрим на результат.

ggplot(data, aes(x = education, y = income, color = type)) +
  geom_point() +
  geom_abline(intercept = coeff[2], slope = coeff[1], color='darkred') +
  geom_abline(intercept = coeff[4], slope = coeff[1], color='darkblue') +
  geom_abline(intercept = coeff[3], slope = coeff[1], color='darkgreen')