Оценка доверительных интервалов 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))

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