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)
## Rows: 300
## Columns: 3
## $ news        <dbl> 25, 13, 21, 18, 10, 11, 29, 51, 23, 11, 22, 35, 77, 37, 4…
## $ old         <dbl> 103, 61, 62, 180, 55, 68, 120, 45, 61, 116, 68, 83, 149, …
## $ mediana_pop <dbl> 5.966097, 5.745970, 5.195969, 5.595928, 5.900597, 6.11865…

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?

Manualmente

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

funcao_theta = function(df) {
  df %>%
    pull(proporcao) %>%
    mean() ## mudar aqui a estatística
}

theta = funcao_theta(lastfm_ajustada)

bootstrap_funcao <- function(x){
  proporcao = x %>% pull(proporcao)
  ## gera uma amostra aleatoria de tamanho da quantidade de linhas dos dados, com reposição, da métrica proporção
  amostra <- sample(proporcao,
                   size = NROW(proporcao),
                   replace = TRUE)
  return(mean(amostra))

}
### theta_chapeu é o estimador do parametro theta (ou seja, theta_chapeu = estimador para média da proporção)
reamostragens = tibble(amostra = 1:4000) %>% # faremos 4000 vezes
  mutate(theta_c = map_dbl(amostra, ~ bootstrap_funcao(lastfm_ajustada)))

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

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

kable(intervalo) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
erro_i erro_s valor_i valor_s
-0.0099264 0.0099499 0.2324396 0.2523158

Como usamos 95% de confiança, essa é a probabilidade da média amostral da proporção estar a uma distância de, no máximo \(\epsilon\) (erro, \(\epsilon\) = 0.0099499), da média populacional. Dessa forma, temos que o intervalo de confiança é dado por [0.2324396, 0.2523158], isso significa que, com 95% de confiança, a nossa média populacional da proporção de novos artistas está dentro desse intervalo. O nosso \(\hat{\theta}\) é 0.242366

Usando a biblioteca

theta <- function(df,i) {
    mean(
        (df %>%
        slice(i) %>%
        pull(proporcao)
    ))
}

## statistic -> métrica que será usada
myBootstrap <- boot(data = lastfm_ajustada, statistic = theta, R=4000)

ci = tidy(myBootstrap, 
          conf.level = .95,
          conf.method = "bca",
          conf.int = TRUE)
glimpse(ci)
## Rows: 1
## Columns: 5
## $ statistic <dbl> 0.242366
## $ bias      <dbl> 5.926423e-05
## $ std.error <dbl> 0.006078198
## $ conf.low  <dbl> 0.2300858
## $ conf.high <dbl> 0.2542181

Usando a biblioteca para calcular o intervalo, obtivemos [0.2310567, 0.2540742] e o \(\theta^\) = 0.242366, exatamente como o calculado manualmente.

Visualizando…

theta = funcao_theta(lastfm_ajustada)

ggplot() +
  geom_rect(
    data = intervalo,
    aes(xmin = valor_i, xmax = valor_s),
    ymin = -Inf,
    ymax = Inf,
    fill = "gold",
    alpha = 1
  ) +
  geom_histogram(
    data = reamostragens,
    aes(theta_c),
   binwidth = .002,
    fill = "white",
    colour = "darkgrey"
  ) + geom_vline(xintercept = theta, color = "blue", size = 0.7) +
  labs(title = expression("Intervalo estimado via bootstrap - Manualmente"))

A variável proporção é contínua, contudo possui muitos valores distintos, então podemos agrupar os dados por classes (intervalos). A base de cada barra representa uma classe e sua altura a frequencia absoluta com que o valor da classe aparece. A partir dessas informações como os valores se espalham ou se distribuem em torno da média e parece que a distribuições dos valores ocorre de maneira simétrica, e nos dá uma pista sobre a natureza da distribuição. O que ocorre é que a amostra terá uma distribuição parecida com a da população de onde foi retirada.

Correlação entre Pop. Mediana e Proporção de artistas novos.

  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 escutado e a proporção dos artistas escutados que eram novos.
media_pop_5 <- lastfm_ajustada %>% filter(mediana_pop > 5 )

Manualmente

funcao_theta_2 = function(df) {
    x <-df %>% mutate(proporcao = news/(news+old)) %>% 
    summarise(correlacao = cor(mediana_pop, proporcao, method= 'pearson')) %>% 
    pull(correlacao)
    return(x)
}
    
theta_pop = funcao_theta_2(media_pop_5)
theta_pop
## [1] -0.05679804
bootstrap_funcao_2 <- function(x){
    boot_x <- sample_n(x, size = NROW(x), replace = TRUE)
    return(funcao_theta_2((boot_x)))
}
    
reamostragens = tibble(amostra = 1:4000) %>% # faremos 4000 vezes
  mutate(theta_c_pop = map_dbl(amostra, ~ bootstrap_funcao_2(media_pop_5)))
confianca = .95
alpha = 1 - confianca
    
intervalo = reamostragens %>%
mutate(erro = theta_c_pop - theta_pop) %>%
summarise(erro_i = quantile(erro, .05),
erro_s = quantile(erro, .95)) %>%
mutate(valor_i = theta_pop + erro_i,
valor_s = theta_pop + erro_s)

intervalo

Como usamos 95% de confiança, essa é a probabilidade da média amostral da proporção estar a uma distância de, no máximo \(\epsilon\) (erro, \(\epsilon\) = 0.1077206), da média populacional. Dessa forma, temos que o intervalo de confiança é dado por [-0.1572458 , 0.05092256], isso significa que, com 95% de confiança, a nossa média populacional da correlação entre popularidade mediana dos artistas escutado e a proporção dos artistas escutados que eram novos está dentro desse intervalo. O nosso \(\hat{\theta}\) é -0.05679804

Usando a biblioteca

theta_pop <- function(df,i) {
     df = df %>%
        filter(mediana_pop > 5) %>%
        slice(i) %>%
        mutate(prop = news/(news + old),
               cor = cor(mediana_pop, prop))
    mean(df$cor)
}
booted_pop <- boot(data = lastfm, 
               statistic = theta_pop, 
               R = 4000)
ci = tidy(booted_pop, 
          conf.level = .95,
          conf.method = "bca",
          conf.int = TRUE)
glimpse(ci)
## Rows: 1
## Columns: 5
## $ statistic <dbl> -0.05679804
## $ bias      <dbl> 0.001705432
## $ std.error <dbl> 0.06375335
## $ conf.low  <dbl> -0.1837948
## $ conf.high <dbl> 0.06265423

Já usando a biblioteca, chegamos a um resultado um pouco diferente. O intervalo de confiança nesse caso é dado por [ -0.1837948 , 0.06265423]. Contudo o chegamos ao mesmo resultado para nosso \(\hat{\theta}\) é -0.05679804

ci %>%
    ggplot(aes(
        x = "",
        y = statistic,
        ymin = conf.low,
        ymax = conf.high
    )) +
    geom_pointrange() +
    geom_point(size = 3) + 
    labs(x = "Correlação entre popularidade mediana dos artistas x proporção de descobertas", 
         y = "")

Visualizando

IC

A correlação entre a popularidade mediana dos artistas e a proporção de novos artistas é muito baixa, ficando entre -0.1 e 0.0.