Os dados

set.seed(12345)

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

lastfm = lastfm %>% 
  sample_n(300) %>% 
  select(news, old, mediana_pop)

glimpse(lastfm)
## Observations: 300
## Variables: 3
## $ news        <dbl> 28, 35, 13, 24, 14, 17, 13, 21, 34, 55, 10, 33, 10, 217...
## $ old         <dbl> 61, 194, 70, 96, 130, 67, 106, 123, 76, 78, 76, 116, 11...
## $ mediana_pop <dbl> 6.105585, 5.376812, 5.713082, 4.564335, 5.782320, 5.532...

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 em geral escutada por usuários?

Para saber a proporção de novos artistas é só calcular o número de artistas novos (news) dividido pelo total de artistas (old + news).

Bootstrap Manual

theta_prop <- function(data){
  x = data %>% mutate(prop = news/(news + old)) %>% pull(prop)
  mean(x)
}

theta_mean <- theta_prop(lastfm)
theta_mean
## [1] 0.2483568
repeticoes = 10000

bootstrap_manual <- function(data){
  prop = data %>% mutate(prop = news/(old + news)) %>% pull(prop)
  boot_x <- sample(prop,              # amostra
                   size = NROW(prop), # tamanho igual ao da entrada
                   replace = TRUE)    # bootstrap
  
  return(mean(boot_x))
}


reamostragem = tibble(i = 1:repeticoes) %>% 
  mutate(bs_mean = map_dbl(i, ~ bootstrap_manual(lastfm)))

reamostragem
interval = reamostragem %>% mutate(err = bs_mean - theta_mean) %>% 
  summarise(err_inf = quantile(err, .05), 
            err_sup = quantile(err, .95))

interval = interval %>% 
  mutate(lim_inf = theta_mean + err_inf, 
         lim_sup = theta_mean + err_sup)

interval

Bootstrap com biblioteca

library(boot)
library(broom)

theta_m <- function(data, n) {
  x = data %>% slice(n) %>% mutate(prop = news/(news + old)) %>% pull(prop)
  mean(x)
}

booted <- boot(data = lastfm,
               statistic = theta_m,
               R = repeticoes)

ic = tidy(booted, 
          conf.level = .95,
          conf.method = "bca",
          conf.int = TRUE)

glimpse(ic)
## Observations: 1
## Variables: 5
## $ statistic <dbl> 0.2483568
## $ bias      <dbl> 6.205002e-06
## $ std.error <dbl> 0.006706292
## $ conf.low  <dbl> 0.2358745
## $ conf.high <dbl> 0.2619652
ggplot() +
    geom_rect(data = interval, aes(xmin = lim_inf, xmax = lim_sup),
    ymin = -Inf, ymax = Inf, fill = "red", alpha = .5) +
  geom_histogram(data = reamostragem,
    aes(bs_mean), binwidth = .005, fill = "grey", colour = "black") +
  geom_vline(xintercept = theta_mean, color = "blue", size = 0.8) +
  labs(y = "Densidade", x = "Média", title = "IC - Bootstrap manual")

ggplot() +
    geom_rect(data = ic, aes(xmin = conf.low, xmax = conf.high),
    ymin = -Inf, ymax = Inf, fill = "red", alpha = .5) +
  geom_histogram(data = reamostragem,
    aes(bs_mean), binwidth = .005, fill = "grey", colour = "black") +
  geom_vline(xintercept = ic$statistic, color = "blue", size = 0.8) +
  labs(y = "Densidade", x = "Média", title = "IC - Bootstrap com biblioteca")

A média (theta) com o uso dos dois bootstrap’s foi exatamente igual 0.2483568. O intervalo de confiança (IC) com o bootstrap manual foi [0.2375502, 0.2596061] e com o bootstrap da biblioteca foi [0.2358228, 0.2621284]. Com 95% de confiança é possível afirmar 24.83% dos novos artistas são geralmente escutados no Last.fm por usuários, com uma margem de erro que varia pouco mais de 1% para o bootstrap manual e menos de 1% para o bootstrap da biblioteca.

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

Crie intervalos com 95% de confiança.

corr = lastfm %>% filter(mediana_pop > 5) %>%
                  mutate(prop = news/(news + old)) %>%
    summarise(pearson = cor(mediana_pop, prop, method = "pearson"),
              kendall = cor(mediana_pop, prop, method = "kendall"),
              spearman = cor(mediana_pop, prop, method = "spearman"))
corr

Podemos ver que há uma correlação fraca negativa entre as variáveis.

Bootstrap manual

theta_corr <- function(data){
  x = data %>% filter(mediana_pop > 5) %>%
      mutate(prop = news/(news + old)) %>%
      summarise(corr = cor(mediana_pop, prop, method = 'pearson')) %>% 
      pull(corr)
  return(x)
}

theta_cor <- theta_corr(lastfm)
theta_cor
## [1] -0.088961
bootstrap_manual_corr <- function(data){
  boot_x <- sample_n(data,              # amostra
                   size = NROW(data), # tamanho igual ao da entrada
                   replace = TRUE)    # bootstrap
  
  return(theta_corr(boot_x))
}


reamostragem_corr = tibble(i = 1:repeticoes) %>% 
  mutate(bs_mean_corr = map_dbl(i, ~ bootstrap_manual_corr(lastfm)))

reamostragem_corr
interval_corr = reamostragem_corr %>% mutate(err = bs_mean_corr - theta_cor) %>% 
  summarise(err_inf = quantile(err, .05), 
            err_sup = quantile(err, .95))

interval_corr = interval_corr %>% 
  mutate(lim_inf = theta_cor + err_inf, 
         lim_sup = theta_cor + err_sup)

interval_corr
theta_c <- function(data, n) {
  data = data %>% slice(n) %>% filter(mediana_pop > 5) %>% mutate(prop = news/(news + old))
  cor(data$mediana_pop, data$prop, method = 'pearson')
}

booted <- boot(data = lastfm,
               statistic = theta_c,
               R = repeticoes)

ic_corr = tidy(booted, 
          conf.level = .95,
          conf.method = "bca",
          conf.int = TRUE)

glimpse(ic_corr)
## Observations: 1
## Variables: 5
## $ statistic <dbl> -0.088961
## $ bias      <dbl> 0.003188351
## $ std.error <dbl> 0.06919008
## $ conf.low  <dbl> -0.2304036
## $ conf.high <dbl> 0.04112106
ggplot() +
    geom_rect(data = interval_corr, aes(xmin = lim_inf, xmax = lim_sup),
    ymin = -Inf, ymax = Inf, fill = "red", alpha = .5) +
  geom_histogram(data = reamostragem_corr,
    aes(bs_mean_corr), binwidth = .005, fill = "grey", colour = "black") +
  geom_vline(xintercept = theta_cor, color = "blue", size = 0.8) +
  labs(y = "Densidade", x = "Correlação", title = "IC - Bootstrap manual")

ggplot() +
    geom_rect(data = ic_corr, aes(xmin = conf.low, xmax = conf.high),
    ymin = -Inf, ymax = Inf, fill = "red", alpha = .5) +
  geom_histogram(data = reamostragem_corr,
    aes(bs_mean_corr), binwidth = .005, fill = "grey", colour = "black") +
  geom_vline(xintercept = ic_corr$statistic, color = "blue", size = 0.8) +
  labs(y = "Densidade", x = "Correlação", title = "IC - Bootstrap com biblioteca")

A correlação (theta_cor) com o uso dos dois bootstrap’s foi exatamente igual -0.088961. O intervalo de confiança (IC) com o bootstrap manual foi [-0.2001402, 0.02894353] e com o bootstrap da biblioteca foi [-0.2291493, 0.04672072]. Com 95% de confiança é possível afirmar que a correlação entre a popularidade mediana dos artistas escutados e a proporção dos artistas escutados que eram novos é de -8.89%, com uma margem de erro que varia de 11% para o bootstrap manual e de 6% para o bootstrap da biblioteca, o que indica que dificilmente essas variáveis estão relacionadas de alguma forma.