Оценка доверительных интервалов bootstrap на примере суперкубка #TidyTuesday

Данная заметка - это любительский перевод статьи Julia Silge.

Эта статья взята из блога Julia Silge, которая демонстрирует использование пакетов tidymodels. В сегодняшней заметке будет показана относительно новая функция из пакета rsample - reg_intervals. Данная функция разработана для быстрого поиска доверительных интервалов bootstrap.

Данные: набор #TidyTuesday о рекламных роликах суперкубка.

Знакомимся с данными

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

Прочитаем наши данные:

library(tidyverse)
youtube <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2021/2021-03-02/youtube.csv")

Давайте сделаем некоторые преобразование и посмотрим, как ведут себя данные с течением времени.

youtube %>%
  select(year, funny:use_sex) %>% 
  pivot_longer(funny:use_sex) %>% 
  group_by(year, name) %>%
  summarise(prop = mean(value)) %>% 
  ungroup() %>%
  ggplot(aes(year, prop, color = name)) +
  geom_line(size = 1.2, show.legend = FALSE) +
  facet_wrap(~name) +
  scale_y_continuous(labels = scales::percent) +
  labs(x = NULL, y = "% of commercials")

Сделаем простую модель

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

simple_mod <- lm(year ~ funny + show_product_quickly + patriotic + celebrity + danger + animals + use_sex,
                 data = youtube)
summary(simple_mod)

Call:
lm(formula = year ~ funny + show_product_quickly + patriotic + 
    celebrity + danger + animals + use_sex, data = youtube)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.5254  -4.1023   0.1456   3.9662  10.1727 

Coefficients:
                          Estimate Std. Error  t value Pr(>|t|)    
(Intercept)              2011.0838     0.9312 2159.748  < 2e-16 ***
funnyTRUE                  -2.8979     0.8593   -3.372  0.00087 ***
show_product_quicklyTRUE    0.7706     0.7443    1.035  0.30160    
patrioticTRUE               2.0455     1.0140    2.017  0.04480 *  
celebrityTRUE               2.4416     0.7767    3.144  0.00188 ** 
dangerTRUE                  0.4814     0.7846    0.614  0.54007    
animalsTRUE                 0.1082     0.7330    0.148  0.88274    
use_sexTRUE                -2.4041     0.8175   -2.941  0.00359 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.391 on 239 degrees of freedom
Multiple R-squared:  0.178,	Adjusted R-squared:  0.1539 
F-statistic: 7.393 on 7 and 239 DF,  p-value: 4.824e-08

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

library(rsample)
bootstraps(youtube, times = 1000)

# Bootstrap sampling 
# A tibble: 1,000 × 2
   splits           id           
   <list>           <chr>        
 1 <split [247/87]> Bootstrap0001
 2 <split [247/85]> Bootstrap0002
 3 <split [247/88]> Bootstrap0003
 4 <split [247/83]> Bootstrap0004
 5 <split [247/87]> Bootstrap0005
 6 <split [247/84]> Bootstrap0006
 7 <split [247/92]> Bootstrap0007
 8 <split [247/84]> Bootstrap0008
 9 <split [247/91]> Bootstrap0009
10 <split [247/86]> Bootstrap0010
# … with 990 more rows

Нам бы пришлось выполнить достаточно много предварительных шагов, чтобы найти доверительные интервалы. Поэтому в пакете rsample предусмотрена функция reg_intervals(), которая находит доверительные интервалы для моделей типа: lm(), glm().

set.seed(123)
youtube_intervals <- reg_intervals(year ~ funny + show_product_quickly + patriotic + celebrity + danger + animals + use_sex,
                                   data = youtube,
                                   type = "percentile",
                                   keep_reps = TRUE)
youtube_intervals

# A tibble: 7 × 7
  term                     .lower .estimate .upper .alpha .method      .replicates
  <chr>                     <dbl>     <dbl>  <dbl>  <dbl> <chr>      <list<tibble>
1 animalsTRUE              -1.22      0.144  1.51    0.05 percentile   [2,001 × 2]
2 celebrityTRUE             0.828     2.46   4.06    0.05 percentile   [2,001 × 2]
3 dangerTRUE               -1.01      0.515  2.09    0.05 percentile   [2,001 × 2]
4 funnyTRUE                -4.58     -2.91  -1.26    0.05 percentile   [2,001 × 2]
5 patrioticTRUE             0.112     2.05   3.88    0.05 percentile   [2,001 × 2]
6 show_product_quicklyTRUE -0.839     0.740  2.23    0.05 percentile   [2,001 × 2]
7 use_sexTRUE              -4.04     -2.43  -0.952   0.05 percentile   [2,001 × 2]

Посмотри на результаты bootstrap

Мы можем визуализировать наши результаты. Обратим внимание, что если бы мы не установил параметр keep_reps = TRUE, то у нас были бы только интервалы и мы бы не смогли построить график №3.

youtube_intervals %>%
  mutate(term = str_remove(term, "TRUE"),
         term = fct_reorder(term, .estimate)) %>%
  ggplot(aes(.estimate, term)) +
  geom_vline(xintercept = 0, size = 1.5, lty = 2, color = "gray80") +
  geom_errorbarh(aes(xmin = .lower, xmax = .upper),
                 size = 1.5, alpha = 0.5, color = "midnightblue") +
  geom_point(size = 3, color = "midnightblue") +
  labs(x = "Increase in year for each commercial characteristic",
       y = NULL)

outube_intervals %>%
  mutate(term = str_remove(term, "TRUE"),
         term = fct_reorder(term, .estimate)) %>%
  unnest(.replicates) %>% 
  ggplot(aes(estimate, fill = term)) +
  geom_vline(xintercept = 0, size = 1.5, lty = 2, color = "gray50") +
  geom_histogram(alpha = 0.8, show.legend = FALSE) +
  facet_wrap(vars(term))

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