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)

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.

Proporção de Novos Artistas

1. Qual a proporção de novos artistas em geral escutada por usuários?

Para calcular esta métrica, fazemos:

lastfm <- lastfm %>%
mutate(prop_new = news / (news + old))

Realizando o processo de bootstrap, teremos as distribuições:

Bootstrap

theta_c <- lastfm %>%
pull(prop_new) %>%
mean()

repeticoes = 4000 # pelo menos 2000, mas mais não faz mal.

bootstrap_mean <- function(x){
  prop_new = x %>% pull(prop_new)
  resample <- sample(prop_new, # amostre dos dados
  size = NROW(prop_new), # tamanho igual ao recebido
  replace = TRUE) # aqui é o bootstrap
  
  return(mean(resample))
}

set.seed(1212)

# A REAMOSTRAGEM
reamostragens = tibble(i = 1:repeticoes) %>%
mutate(theta_c_s = map_dbl(i, ~ bootstrap_mean(lastfm)))

Distribuições das médias amostrais

reamostragens %>%
  ggplot(aes(x = theta_c_s)) +
  geom_histogram(colour = "darkorange",
  fill = "white") +
  labs(title = "Distrbuição amostral")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

reamostragens %>%
  ggplot(aes(x = theta_c_s - theta_c)) +
  geom_histogram(colour = "darkblue",
  fill = "white") + 
  labs(title = "Distribuição do erro amostral")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

A partir destas, podemos calcular os intervalos de confiança.

ICs

intervalo = reamostragens %>%
  mutate(erro = theta_c_s - theta_c) %>%
  summarise(erro_i = quantile(erro, .025),
            erro_s = quantile(erro, .975))

intervalo = intervalo %>%
  mutate(valor_i = theta_c + erro_i,
         valor_s = theta_c + erro_s)

ggplot() +
  geom_rect(data = intervalo, aes(xmin = valor_i, xmax = valor_s), 
            ymin = -Inf, ymax = Inf, fill = "gold", alpha = .25) +
  geom_histogram(data = reamostragens, aes(theta_c_s), fill = "white",
                 colour = "darkgrey") +
  geom_vline(xintercept = theta_c, color = "dark green") +
  labs(title = expression("Intervalo estimado via bootstrap"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Por fim, em 95% das reamostragens, a proporção de artistas novos escutada pelos usuários esteve no intervalo [0.2308375, 0.2544643].

ICs usando a biblioteca Boot

Iremos agora realizar o bootstrap com o auxílio de uma biblioteca chamada Boot

library(boot)
library(broom)

bootstrap_mean <- function(data, indices){
  d <- data[indices,]
  
  prop_new = d %>% pull(prop_new)
  
  return(mean(prop_new))
}

booted_mean <- boot(data = lastfm, 
               statistic = bootstrap_mean, 
               R = 2000)

ci = tidy(booted_mean, 
          conf.level = .95,
          conf.method = "bca",
          conf.int = TRUE)

ggplot() +
  geom_rect(data = intervalo, aes(xmin = valor_i, xmax = valor_s), 
            ymin = -Inf, ymax = Inf, fill = "gold", alpha = .25) +
  geom_rect(aes(xmin = ci$conf.low, xmax = ci$conf.high),
            ymin = -Inf, ymax = Inf, fill = "blue", alpha = .25) +
  geom_histogram(data = reamostragens, aes(theta_c_s), fill = "white",
                 colour = "darkgrey") +
  geom_vline(xintercept = theta_c, color = "dark green") +
  labs(title = expression("Intervalo estimado via bootstrap"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Os valores do IC obtidos através da biblioteca Boot são portanto muito próximos aos nossos.

IC calculado: [0.2308375, 0.2544643]

IC usando Boot: [0.2305321, 0.2543645]

Correlação entre popularidade e proporção

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 escutados e a proporção dos artistas escutados que eram novos.

lastfm_pop <- lastfm %>%
  filter(mediana_pop > 5)

Bootstrap

prop_new <- lastfm_pop %>% 
  pull(prop_new)

mediana_pop <- lastfm_pop %>%
  pull(mediana_pop)

theta_c_cor = cor(prop_new, mediana_pop)

repeticoes = 4000 # pelo menos 2000, mas mais não faz mal.

um_bootstrap_cor <- function(x){
  
  resample <- sample_n(x, # amostre dos dados
                     size = NROW(prop_new), # tamanho igual ao recebido
                     replace = TRUE) # aqui é o bootstrap
  
  prop_new = resample %>% pull(prop_new)
  mediana_pop = resample %>% pull(mediana_pop)
                      
  return(cor(prop_new, mediana_pop))
}

set.seed(1212)

# A REAMOSTRAGEM
reamostragens_cor = tibble(i = 1:repeticoes) %>%
mutate(theta_c_s = map_dbl(i, ~ um_bootstrap_cor(lastfm_pop)))

Distribuições das médias amostrais

reamostragens_cor %>%
  ggplot(aes(x = theta_c_s)) +
  geom_histogram(colour = "darkorange",
  fill = "white") +
  labs(title = "Distrbuição amostral")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

  reamostragens_cor %>%
  ggplot(aes(x = theta_c_s - theta_c)) +
  geom_histogram(colour = "darkblue",
  fill = "white") + 
  labs(title = "Distribuição do erro amostral")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

A distribuição amostral a correlação entre as duas variáveis, parece sugerir uma correlação negativa. Vamos calcular os intervalos de confiança para ter maior certeza.

ICs

intervalo_cor = reamostragens_cor %>%
  mutate(erro = theta_c_s - theta_c) %>%
  summarise(erro_i = quantile(erro, .025),
            erro_s = quantile(erro, .975))

intervalo_cor = intervalo_cor %>%
  mutate(valor_i = theta_c + erro_i,
         valor_s = theta_c + erro_s)

ggplot() +
  geom_rect(data = intervalo_cor, aes(xmin = valor_i, xmax = valor_s), 
            ymin = -Inf, ymax = Inf, fill = "gold", alpha = .25) +
  geom_histogram(data = reamostragens_cor, aes(theta_c_s), fill = "white",
                 colour = "darkgrey") +
  geom_vline(xintercept = theta_c_cor, color = "dark green") +
  labs(title = expression("Intervalo estimado via bootstrap"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Em 95% das reamostragens, a correlação entre a popularidade mediana dos artistas escutado e a proporção dos artistas escutados que eram novos esteve no intervalo [-0.1799214, 0.0678431].

Apesar da maior parte do intervalo de confiança estar contida em um valor abaixo de zero, o intervalo contém não só o zero mas também alguns valores positivos. Desta forma não podemos afirmar com confiança, que a correlação seja negativa, apenas que é fraca. pois em 95% das reamostragens, tivemos correlações que não atingiram valores significativos, próximos a 1 ou -1.

ICs usando a biblioteca Boot

bootstrap_cor <- function(data, indices){
  d <- data[indices,]
  
  prop_new = d %>% pull(prop_new)
  mediana_pop = d %>% pull(mediana_pop)
  
  return(cor(prop_new, mediana_pop))
}

booted_cor <- boot(data = lastfm_pop, 
               statistic = bootstrap_cor, 
               R = 2000)

ci_cor = tidy(booted_cor, 
          conf.level = .95,
          conf.method = "bca",
          conf.int = TRUE)

ggplot() +
  geom_rect(data = intervalo_cor, aes(xmin = valor_i, xmax = valor_s), 
            ymin = -Inf, ymax = Inf, fill = "gold", alpha = .25) +
  geom_rect(aes(xmin = ci_cor$conf.low, xmax = ci_cor$conf.high),
            ymin = -Inf, ymax = Inf, fill = "blue", alpha = .25) +
  geom_histogram(data = reamostragens_cor, aes(theta_c_s), fill = "white",
                 colour = "darkgrey") +
  geom_vline(xintercept = theta_c_cor, color = "dark green") +
  labs(title = expression("Intervalo estimado via bootstrap"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Mais uma vez temos valores bem parecidos usando as duas formas de calcular os intervalos.

IC calculado: [-0.1799214, 0.0678431]

IC usando Boot: [-0.1759809, 0.0587894]