Заметки по R: Апельсины и простая регрессия

Подгрузим необходимые пакеты:

library("tidyverse") # коллекция пакетов: ggplot2, dplyr, ...
library("lmtest") # тесты для линейных моделей
library("memisc") # сравнение моделей
library("pander") # таблички в markdown
library("broom") # стандартизация информации о модели 
library("psych") # описательные статистики
library("modelr") # добавление прогнозов/остатков

Простая работа с данными до оценки регрессий

В R есть куча встроенных наборов данных. Многие пакеты в себе помимо функций содержат также ещё и примеры данных. Например, пакет ggplot2 содержит в себе данные по стоимости бриллиантов diamonds.

Список всех наборов данных можно узнать с помощью команды data().

Мы воспользуемся набором данных Orange про апельсиновые деревья.

Совсем кратко о наборе данных:

glance(Orange)

##   nrow ncol complete.obs na.fraction
## 1   35    3           35           0

Описательные статистики:

tidy(Orange) 

##          column  n     mean        sd median  trimmed      mad min  max
## 1         Tree* 35   3.0000   1.43486      3   3.0000   1.4826   1    5
## 2           age 35 922.1429 491.86453   1004 937.0690 545.5968 118 1582
## 3 circumference 35 115.8571  57.48818    115 115.1379  77.0952  30  214
##   range          skew  kurtosis         se
## 1     4  0.000000e+00 -1.395755  0.2425356
## 2  1464 -2.556469e-01 -1.293623 83.1402798
## 3   184  1.655431e-05 -1.239487  9.7172759

Как выглядит начало таблички?

head(Orange)

## Grouped Data: circumference ~ age | Tree
##   Tree  age circumference
## 1    1  118            30
## 2    1  484            58
## 3    1  664            87
## 4    1 1004           115
## 5    1 1231           120
## 6    1 1372           142

Структура набора данных.

str(Orange)

## Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame':   35 obs. of  3 variables:
##  $ Tree         : Ord.factor w/ 5 levels "3"<"1"<"5"<"2"<..: 2 2 2 2 2 2 2 4 4 4 ...
##  $ age          : num  118 484 664 1004 1231 ...
##  $ circumference: num  30 58 87 115 120 142 145 33 69 111 ...
##  - attr(*, "formula")=Class 'formula'  language circumference ~ age | Tree
##   .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv> 
##  - attr(*, "labels")=List of 2
##   ..$ x: chr "Time since December 31, 1968"
##   ..$ y: chr "Trunk circumference"
##  - attr(*, "units")=List of 2
##   ..$ x: chr "(days)"
##   ..$ y: chr "(mm)"

Команда str — рабочая лошадка! Пакеты R пишет огромное количество людей, а она позволяет понять, из чего состоит тот или иной объект.

Более подробную информацию о наборе данных про Апельсины ищем у Чебурашки. Чебурашку можно позвать командой help(Orange).

Графики

ggplot(Orange) + 
  geom_point(aes(x = age, y = circumference, color = Tree)) + 
  labs(x = "Возраст дерева, лет", 
       y = "Диаметр ствола, мм",
  title = "Пять апельсиновых деревьев") +
  scale_colour_discrete(name = "Дерево п/п")

Оценим обыкновенную парную регрессию для того, чтобы понять, как диаметр дерева зависит от его возраста:

model <- lm(circumference ~ age, data = Orange) # оцениваем модель
report <- summary(model) # создаём отчет по модели
report # выводим отчет на экран

## 
## Call:
## lm(formula = circumference ~ age, data = Orange)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -46.310 -14.946  -0.076  19.697  45.111 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 17.399650   8.622660   2.018   0.0518 .  
## age          0.106770   0.008277  12.900 1.93e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.74 on 33 degrees of freedom
## Multiple R-squared:  0.8345, Adjusted R-squared:  0.8295 
## F-statistic: 166.4 on 1 and 33 DF,  p-value: 1.931e-14

Из результатов оценивания можно вытащить

Ковариационную матрицу, \(\hat{\sigma}^2(X'X)^{-1}\):

vcov(model)

##             (Intercept)           age
## (Intercept) 74.35026205 -6.316908e-02
## age         -0.06316908  6.850249e-05

Любимый всеми \(R^2\):

report$r.squared

## [1] 0.8345167

А также скорректированный \(R^2_{adj}\):

report$adj.r.squared

## [1] 0.829502

Коэффициенты:

coef(model)

## (Intercept)         age 
##  17.3996502   0.1067703

Доверительные интервалы для коэффициентов:

confint(model, level = 0.90)

##                   5 %       95 %
## (Intercept) 2.8070030 31.9922974
## age         0.0927633  0.1207774

Отдельно табличка с результатами тестов:

coeftest(model)

## 
## t test of coefficients:
## 
##               Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) 17.3996502  8.6226598  2.0179   0.05179 .  
## age          0.1067703  0.0082766 12.9002 1.931e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Несколько моделей можно сравнить рядом в одной табличке:

model_2 <- lm(data = Orange, 
              circumference ~ age + I(age^2))
mtable(model, model_2)

## 
## Calls:
## model: lm(formula = circumference ~ age, data = Orange)
## model_2: lm(formula = circumference ~ age + I(age^2), data = Orange)
## 
## =========================================
##                     model     model_2    
## -----------------------------------------
##   (Intercept)     17.400      10.287     
##                   (8.623)    (13.199)    
##   age              0.107***    0.131***  
##                   (0.008)     (0.035)    
##   I(age^2)                    -0.000     
##                               (0.000)    
## -----------------------------------------
##   R-squared            0.8        0.8    
##   adj. R-squared       0.8        0.8    
##   sigma               23.7       23.9    
##   F                  166.4       82.2    
##   p                    0.0        0.0    
##   Log-likelihood    -159.5     -159.2    
##   Deviance         18594.7    18301.8    
##   AIC                325.0      326.4    
##   BIC                329.6      332.6    
##   N                   35         35      
## =========================================

Построим точечный прогноз диаметра дерева:

new.data <- data_frame(age = c(100, 200, 300))
predict(model, new.data)

##        1        2        3 
## 28.07668 38.75372 49.43075

Прогноз с доверительным интервалом для среднего диаметра дерева:

predict(model, new.data, interval = "confidence")

##        fit      lwr      upr
## 1 28.07668 12.00511 44.14826
## 2 38.75372 24.10765 53.39979
## 3 49.43075 36.14955 62.71195

И прогноз с предиктивным интервалом для конкретного наблюдения:

predict(model, new.data, interval = "prediction")

##        fit         lwr      upr
## 1 28.07668 -22.8219365 78.97530
## 2 38.75372 -11.7129209 89.22035
## 3 49.43075  -0.6568182 99.51831

Есть удобный пакет broom. Основная его философия проста — всю инфу о модели можно выдавать в data.frame. Пакет broom делит информацию на три уровня:

• краткая о модели в целом
• информация о коэффициентах
• то, что можно добавить к исходным данным после оценки модели.

О модели в целом:

glance(model)

##   r.squared adj.r.squared    sigma statistic      p.value df    logLik
## 1 0.8345167      0.829502 23.73767  166.4159 1.930596e-14  2 -159.4804
##        AIC      BIC deviance df.residual
## 1 324.9607 329.6268 18594.74          33

Более подробную о коэффициентах:

tidy(model)

##          term   estimate   std.error statistic      p.value
## 1 (Intercept) 17.3996502 8.622659801  2.017898 5.179267e-02
## 2         age  0.1067703 0.008276623 12.900228 1.930596e-14

Довесим исходные данные предсказанными значениями, остатками и прочим

orange_plus <- augment(model, Orange)
glimpse(orange_plus)

## Observations: 35
## Variables: 10
## $ Tree          <ord> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, ...
## $ age           <dbl> 118, 484, 664, 1004, 1231, 1372, 1582, 118, 484,...
## $ circumference <dbl> 30, 58, 87, 115, 120, 142, 145, 33, 69, 111, 156...
## $ .fitted       <dbl> 29.99855, 69.07649, 88.29515, 124.59706, 148.833...
## $ .se.fit       <dbl> 7.771498, 5.408300, 4.545789, 4.069196, 4.757519...
## $ .resid        <dbl> 0.001451402, -11.076487573, -1.295146086, -9.597...
## $ .hat          <dbl> 0.10718481, 0.05190932, 0.03667265, 0.02938603, ...
## $ .sigma        <dbl> 24.10572, 24.02169, 24.10459, 24.04413, 23.53757...
## $ .cooksd       <dbl> 2.513501e-10, 6.286999e-03, 5.882009e-05, 2.5492...
## $ .std.resid    <dbl> 0.0000647096, -0.4792244734, -0.0555896161, -0.4...

К исходным данным добавлены:

• .fitted — предсказанные \(\hat y_i\)
• .se.fit — стандартное отклонение прогнозов \(se(\hat y_i)\)
• .resid — остатки \(\hat\varepsilon_i\)
• .hat — диагональный элемент матрицы-шляпницы (hat matrix), \(\hat{y}=Hy\), \(H=X(X'X)^{-1}X'\), вес, с которым \(y_i\) входит в формулу для \(\hat y_i\).
• .sigma — какой станет оценка \(\hat\sigma\), если выкинуть \(i\)-ое наблюдение
• .std.resid — стандартизованные остатки, \(s_i = \frac{\hat\varepsilon_i}{\sqrt{\hat{\sigma}^2(1-h_{ii})}}\)
• .cooksd — расстояние Кука, \(D_i\), \[ D_i = \frac{ \sum_{j=1}^n (\hat y_j\ - \hat y_{j(-i)})^2 }{k \ \mathrm{MSE}} \]

Однако чаще всего вместо всего этого барахла нужны лишь прогнозы:

orange_plus_fitted <- Orange %>% add_predictions(model)
glimpse(orange_plus_fitted)

## Observations: 35
## Variables: 4
## $ Tree          <ord> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, ...
## $ age           <dbl> 118, 484, 664, 1004, 1231, 1372, 1582, 118, 484,...
## $ circumference <dbl> 30, 58, 87, 115, 120, 142, 145, 33, 69, 111, 156...
## $ pred          <dbl> 29.99855, 69.07649, 88.29515, 124.59706, 148.833...

Или в крайнем случае остатки:

orange_plus_residuals <- Orange %>% add_residuals(model)
glimpse(orange_plus_residuals)

## Observations: 35
## Variables: 4
## $ Tree          <ord> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, ...
## $ age           <dbl> 118, 484, 664, 1004, 1231, 1372, 1582, 118, 484,...
## $ circumference <dbl> 30, 58, 87, 115, 120, 142, 145, 33, 69, 111, 156...
## $ resid         <dbl> 0.001451402, -11.076487573, -1.295146086, -9.597...

Расстояние Кука

Но почему аборигены съели Кука? За что — неясно, молчит наука. Владимир Высоцкий

Расстояние Кука показывает, насколько сильно данное наблюдение влияет на оценённую регрессию. Насколько сильно изменятся прогнозы и оценки коэффициентов, если мы выкинем наблюдение номер \(i\)? Есть у него два эквивалентных определения: \[ D_i = \frac{ \left(\hat{\beta} - \hat{\beta}^{(-i)} \right)' X'X \left(\hat{\beta} - \hat{{\beta}}^{(-i)}\right) } {(1+k)\hat\sigma^2} \] И: \[ D_i = \frac{\hat \varepsilon_i^2}{k \ \mathrm{MSE}}\left[\frac{h_{ii}}{(1-h_{ii})^2}\right] \] Здесь:

• \(MSE\) — среднеквадратичная ошибка, \(MSE=RSS/n=\sum (y_i -\hat y_i)^2/n\)
• \(\hat y_j\) — прогноз \(j\)-го наблюдения, построенный по всем наблюдениям
• \(\hat y_{j(-i)}\) — прогноз \(j\)-го наблюдения, построенный по всем наблюдениям кроме \(i\)-го
• \(\hat \beta\) — оценки коэффициентов, построенные по всем наблюдениям
• \(\hat \beta^{(-i)}\) — оценки коэффициентов, построенные по всем наблюдениям кроме \(i\)-го

С помощью расстояния Кука можно обнаруживать наблюдения сильно влияющие на результат, выбросы.

orange_plus <- mutate(orange_plus, i = row_number())
ggplot(data = orange_plus) + 
  geom_bar(aes(x = i, 
               y = .cooksd,
               fill = age),
           stat = "identity") + 
  labs(x = "Номер наблюдения",
       y = "Расстояние Кука",
       title = "Влиятельность наблюдений")

На картинке мы видим, что в целом наибольшее влияние на нашу модель оказывают взрослые деревья. Два самах больших расстояния Кука — у очень старых деревьев. Напомним, что наблюдения отсортированы по номеру дерева, а для каждого дерева — по возрасту дерева.

Что ещё можно вытащить из оцененной модели, можно узнать несколькими способами. Например, стоит почитать help(summary.lm). Можно ещё набрать в консоли model$ и пару раз кликнуть <Tab>. ещё имеет смысл попробовать нажать <Tab> после summ$. Можно набрать str(model). Да ещё гугл есть.

А вообще бывает полезно изобрести велосипед. Пусть где-то в R уже есть функция, считающая доверительный интервал для прогноза. Если вместо этого написать эту функцию своими руками, то можно получить плюс один к экспириенсу и перейти на следующий уровень.

Почиташки:

• Виньетка к пакету broom
• Wiki, Cook’s distance
• Herman Aguinis, Best practice for handling outliers