Os dados

lastfm = read_csv(here::here("data/experimento-lastfm.csv"), 
                  col_types = cols(.default = col_double(), 
                                   user = col_character()))

lastfm = lastfm %>% 
  select(news, old, mediana_pop) %>% 
  na.omit(.) %>% 
  mutate(prop_news = news/(news+old))

glimpse(lastfm)
## Observations: 11,988
## Variables: 4
## $ news        <dbl> 38, 32, 24, 36, 28, 25, 30, 28, 14, 22, 27, 21, 23...
## $ old         <dbl> 116, 130, 10, 46, 128, 95, 228, 53, 46, 118, 90, 3...
## $ mediana_pop <dbl> 5.660406, 4.730935, 5.896073, 5.214290, 5.865121, ...
## $ prop_news   <dbl> 0.2467532, 0.1975309, 0.7058824, 0.4390244, 0.1794...

Proporção de artistas novos e popularidade

Utilizaremos ICs para estimar duas métricas sobre os usuários do LastFM em geral durante um período de 6 meses. Em ambos os casos faremos isso a partir de uma amostra de 300 usuários. As duas métricas são:

  1. Qual a proporção de novos artistas escutados pelos usuários
  2. Para os usuários que gostam de música muito pop (mediana_pop > 5), qual a correlação entre a popularidade mediana dos artistas escutado e a proporção dos artistas escutados que eram novos.

Aplicando bootstrapping manual

Proporção de novos artistas

funcao_theta_prop_news = function(df) {
  df %>%
    pull(prop_news) %>%
    mean()
}

theta = funcao_theta_prop_news(lastfm)
summary(lastfm)
##       news             old          mediana_pop      prop_news      
##  Min.   : 10.00   Min.   :   6.0   Min.   :3.327   Min.   :0.01462  
##  1st Qu.: 16.00   1st Qu.:  62.0   1st Qu.:5.360   1st Qu.:0.16239  
##  Median : 24.00   Median :  87.0   Median :5.670   Median :0.22667  
##  Mean   : 32.52   Mean   : 100.3   Mean   :5.579   Mean   :0.24680  
##  3rd Qu.: 40.00   3rd Qu.: 122.0   3rd Qu.:5.884   3rd Qu.:0.31178  
##  Max.   :436.00   Max.   :1213.0   Max.   :6.408   Max.   :0.89831
set.seed(12345)
amostra = lastfm %>%  
  sample_n(300) 

theta_c = funcao_theta_prop_news(amostra)
repeticoes = 5000

um_bootstrap <- function(x){
  prop_news = x %>% pull(prop_news)
  boot_x <- sample(prop_news,
                   size = NROW(x),
                   replace = TRUE)
  return(mean(boot_x))
}

reamostragens = tibble(i = 1:repeticoes) %>% 
  mutate(theta_c_s = map_dbl(i, ~ um_bootstrap(amostra)))
reamostragens %>%
  ggplot(aes(x = theta_c_s)) +
  geom_histogram(binwidth = 0.002,
                 colour = "darkorange",
                 fill = "white")

reamostragens %>%
  ggplot(aes(x = theta_c_s - theta_c)) +
  geom_histogram(binwidth = 0.002,
                 colour = "darkblue",
                 fill = "white")

intervalo = reamostragens %>% 
  mutate(erro = theta_c_s - theta_c) %>% 
  summarise(erro_i = as.numeric(quantile(erro, .05)), 
            erro_s = as.numeric(quantile(erro, .95)))

intervalo
## # A tibble: 1 x 2
##     erro_i  erro_s
##      <dbl>   <dbl>
## 1 -0.00951 0.00998
intervalo = intervalo %>% 
  mutate(valor_i = theta_c + erro_i, 
         valor_s = theta_c + erro_s)

intervalo
## # A tibble: 1 x 4
##     erro_i  erro_s valor_i valor_s
##      <dbl>   <dbl>   <dbl>   <dbl>
## 1 -0.00951 0.00998   0.228   0.248
inf = round(intervalo$valor_i, digits = 3)
sup = round(intervalo$valor_s, digits = 3)
t_c = round(theta_c, digits = 3)
t = round(theta, digits = 3)

ggplot() +
  geom_rect(
    data = intervalo,
    aes(xmin = inf, xmax = sup),
    ymin = -Inf,
    ymax = Inf,
    fill = "gold",
    alpha = .25
  ) +
  geom_histogram(
    data = reamostragens,
    aes(theta_c_s),
    binwidth = .002,
    fill = "white",
    colour = "darkgrey"
  ) +
  geom_vline(xintercept = theta,
             color = "blue",
             size = 1) +
  geom_vline(xintercept = t_c, color = "dark green") +
  scale_x_continuous(breaks = c(inf, sup, t_c)) +
  labs(title = expression("Intervalo estimado via bootstrap"))

Obs: Estamos destacando o theta apenas por curiosidade, considerando que no mundo real não o conhecemos.

A partir desta amostra, estimamos que os usuários escutam, em média, 23,8% novos artistas (95% CI [0.228, 0.248]).

Correlação entre a popularidade mediana dos artistas escutados e a proporção dos artistas escutados que eram novos

funcao_theta_corr = function(df) {
    df <-  df %>% select(prop_news, mediana_pop) 
    t = cor(df$prop_news, df$mediana_pop)
}

theta_corr = funcao_theta_corr(lastfm)
theta_c_corr = funcao_theta_corr(amostra)
repeticoes = 5000

um_bootstrap <- function(x){
  boot_x <- sample_n(x,
                   size = NROW(x), 
                   replace = TRUE) 
  theta_c_s_corr = funcao_theta_corr(boot_x)
  return(theta_c_s_corr)
}

reamostragens = tibble(i = 1:repeticoes) %>% 
  mutate(theta_c_s_corr = map_dbl(i, ~ um_bootstrap(amostra)))
reamostragens %>%
  ggplot(aes(x = theta_c_s_corr)) +
  geom_histogram(binwidth = 0.02,
                 colour = "darkorange",
                 fill = "white")

reamostragens %>%
  ggplot(aes(x = theta_c_s_corr - theta_c)) +
  geom_histogram(binwidth = 0.02,
                 colour = "darkblue",
                 fill = "white")

intervalo = reamostragens %>% 
  mutate(erro = theta_c_s_corr - theta_c_corr) %>% 
  summarise(erro_i = as.numeric(quantile(erro, .05)), 
            erro_s = as.numeric(quantile(erro, .95)))

intervalo
## # A tibble: 1 x 2
##    erro_i erro_s
##     <dbl>  <dbl>
## 1 -0.0889 0.0937
intervalo = intervalo %>% 
  mutate(valor_i = theta_c_corr + erro_i, 
         valor_s = theta_c_corr + erro_s)

intervalo
## # A tibble: 1 x 4
##    erro_i erro_s valor_i valor_s
##     <dbl>  <dbl>   <dbl>   <dbl>
## 1 -0.0889 0.0937 -0.0840  0.0987
inf = round(intervalo$valor_i, digits = 3)
sup = round(intervalo$valor_s, digits = 3)
t_c = round(theta_c_corr, digits = 3)
t = round(theta_corr, digits = 3)

ggplot() +
  geom_rect(
    data = intervalo,
    aes(xmin = inf, xmax = sup),
    ymin = -Inf,
    ymax = Inf,
    fill = "gold",
    alpha = .25
  ) +
  geom_histogram(
    data = reamostragens,
    aes(theta_c_s_corr),
    binwidth = .02,
    fill = "white",
    colour = "darkgrey"
  ) +
  geom_vline(xintercept = t,
             color = "blue",
             size = 1) +
  geom_vline(xintercept = t_c, color = "dark green") +
  scale_x_continuous(breaks = c(inf, sup, t_c, t)) +
  labs(title = expression("Intervalo estimado via bootstrap"))

A partir desta amostra, estimamos uma baixa e positiva correlação entre a popularidade mediana dos artistas escutados e a proporção dos artistas escutados que eram novos (95% CI [ -0.084, 0.099]).

Utilizando a biblioteca

Proporção de novos artistas

set.seed(12345)
func_theta_prop_news = function(df, i) {
    df %>% 
        slice(i) %>% 
        summarise(prop_news = mean(prop_news)) %>% 
        pull(prop_news)
}
amostra = amostra %>% ungroup()
prop_news_boot= amostra %>% boot(statistic = func_theta_prop_news, R = 5000) %>% 
    tidy(conf.level = .95,
         conf.int = TRUE)
prop_news_boot
## # A tibble: 1 x 5
##   statistic      bias std.error conf.low conf.high
##       <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1     0.238 0.0000175   0.00587    0.226     0.249

Temos aqui um resultado muito próximo da abordagem manual.

Correlação entre a popularidade mediana dos artistas escutados e a proporção dos artistas escutados que eram novos

set.seed(12345)
func_theta_corr = function(d, i) {
    d = d %>% 
        slice(i)
    corr = funcao_theta_corr(d)
    return(corr)
}
amostra = amostra %>% ungroup()
corr_boot = amostra %>% boot(statistic = func_theta_corr, R = 5000) %>% 
    tidy(conf.level = .95,
         conf.int = TRUE)
corr_boot
## # A tibble: 1 x 5
##   statistic      bias std.error conf.low conf.high
##       <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1   0.00492 -0.000704    0.0545   -0.101     0.110

Aqui, a diferença cresce um pouco, mas ainda pode ser considerado bem próximo do resultado obtido com a abordagem manual.