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)

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

glimpse(lastfm)
## Observations: 300
## Variables: 4
## $ news        <dbl> 25, 13, 21, 18, 10, 11, 29, 51, 23, 11, 22, 35, 77...
## $ old         <dbl> 103, 61, 62, 180, 55, 68, 120, 45, 61, 116, 68, 83...
## $ mediana_pop <dbl> 5.966097, 5.745970, 5.195969, 5.595928, 5.900597, ...
## $ proporcao   <dbl> 0.19531250, 0.17567568, 0.25301205, 0.09090909, 0....

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?
  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.

Crie intervalos com 95% de confiança.

Primeiramente para responder essas questões, calculamos na amostra utilizada a média da proporção de novos artistas em geral escutada por usuários e a correlação entre a popularidade mediana e a proporção de novos artistas.

funcao_theta = function(df) {
  df %>%
    summarise(mean_prop = mean(proporcao), correlacao = cor(mediana_pop[mediana_pop > 5], proporcao[mediana_pop > 5], method = "pearson" ) )
}

theta = funcao_theta(lastfm)

theta_prop = as.double(theta[1,1])
theta_cor = as.double(theta[1,2])

theta_prop
## [1] 0.242366
theta_cor
## [1] -0.05679804

O segundo passo foi realizar 4000 reamostragens a partir das amostra inicial, e calcular a média da proporção e da correlação.

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

um_bootstrap <- function(df){
  boot_x <- sample_n(df,           # amostre dos dados
                   size = NROW(news), # tamanho igual ao recebido
                   replace = TRUE) # aqui é o bootstrap
  return(funcao_theta(boot_x))
}

um_bootstrap(lastfm)
# A REAMOSTRAGEM
reamostragens = tibble(i = 1:repeticoes) %>% 
  mutate(theta_prop_s = map(i, ~ um_bootstrap(lastfm))) %>%
    unnest()

reamostragens

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

Com os dados anteriormente gerados, temos um histograma da distribuição amostral das médias da proporção das reamostras geradas e um histograma da diferença das médias da proporção das reamostras para a da amostra original.

reamostragens %>%
  ggplot(aes(x = mean_prop)) +
  geom_histogram(binwidth = .002,
                 colour = "darkorange",
                 fill = "white")

reamostragens %>%
  ggplot(aes(x = mean_prop - theta_prop)) +
  geom_histogram(binwidth = .004,
                 colour = "darkorange",
                 fill = "white")

A partir daí, calculamos o erro inferior e o erro superior, levando em conta criar um intervalo com 95% dos valores que a média da proporção pode assumir nas reamostras. Para isso utilizamos o 2.5 percentil e o 97.5 percentil.

intervalo = reamostragens %>% 
  mutate(erro_prop = mean_prop - theta_prop) %>% 
  summarise(erro_prop_i = quantile(erro_prop, .025), 
            erro_prop_s = quantile(erro_prop, .975))
intervalo

E com isso obtemos o intervalo de confiança subtraindo o erro inferior e somando o erro superior da média da amostra original.

intervalo = intervalo %>% 
  mutate(valor_i = theta_prop + erro_prop_i, 
         valor_s = theta_prop + erro_prop_s)

intervalo

Assim temos que a proporção de novos artistas em geral escutada por usuários estimada é 0.242366, ou seja 24.23% dos artistas escutados são novos, podendo variar 0.01153788 para menos, e 0.0118495 para mais. Temos assim um intervalo entre 23.08% e 25.42% , com uma probabilidade de 95% do intervalo refletir a realidade, ou seja, 95% de changes da média real da população está dentro desse intervalo.

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(mean_prop),
    binwidth = .003,
    fill = "white",
    colour = "darkgrey"
  ) +
  geom_vline(xintercept = theta_prop, color = "blue", size = 1.2) +
  labs(title = expression("Intervalo estimado via bootstrap"))

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

Seguimos os passos anteriores para questão 1, e obtemos um histograma da distribuição amostral da correlação entre a popularidade medidana e a proporção de novos artistas escutados das reamostras geradas, e um histograma da diferença dessa estatística das reamostras para a da amostra original.

reamostragens %>%
  ggplot(aes(x = correlacao)) +
  geom_histogram(binwidth = .02,
                 colour = "darkorange",
                 fill = "white")

reamostragens %>%
  ggplot(aes(x = correlacao - theta_cor)) +
  geom_histogram(binwidth = .04,
                 colour = "darkorange",
                 fill = "white")

A partir daí, calculamos o erro inferior e o erro superior, levando em conta criar um intervalo com 95% dos valores que a correlação pode assumir nas reamostras. Para isso utilizamos o 2.5 percentil e o 97.5 percentil.

intervalo_cor = reamostragens %>% 
  mutate(erro_cor = correlacao - theta_cor) %>% 
  summarise(erro_cor_i = quantile(erro_cor, .025), 
            erro_cor_s = quantile(erro_cor, .975))
intervalo_cor

E com isso obtemos o intervalo de confiança.

intervalo_cor = intervalo_cor %>% 
  mutate(valor_i_cor = theta_cor + erro_cor_i, 
         valor_s_cor = theta_cor + erro_cor_s)

intervalo_cor

Assim temos que a correlação (utilizando o método de pearson) entre a popularidade mediana dos artistas escutado e a proporção de novos artista escutados é -0.05679804, podendo variar 0.1240708 para menos, e 0.1239988 para mais. Temos assim um intervalo entre -0.1808688 e 0.06720078 com uma probabilidade de 95% do intervalo refletir a realidade, ou seja, 95% de changes da média real da população está contiga nesse intervalo. Então estimasse que exista uma correlação fraca entre as variáveis que pode ser negativa ou positiva.

ggplot() +
  geom_rect(
    data = intervalo_cor,
    aes(xmin = valor_i_cor, xmax = valor_s_cor),
    ymin = -Inf,
    ymax = Inf,
    fill = "gold",
    alpha = .25
  ) +
  geom_histogram(
    data = reamostragens,
    aes(correlacao),
    binwidth = .03,
    fill = "white",
    colour = "darkgrey"
  ) +
  geom_vline(xintercept = theta_cor, color = "blue", size = 1.2) +
  labs(title = expression("Intervalo estimado via bootstrap"))

Gerando os intervalos de confiança utlizando a biblioteca boot

Utilizando a biblioteca boot para gerar os intervalos de confiança, obtemos resultados bastante parecidos. A diferença é esperada já que a cada execução da reamostragem são geradas reamostras diferentes aleatoriamente.

Para a proporção temos o intervalo [0.2310567, 0.2541075] - [23.10%, 25,41%]

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

booted <- boot(data = lastfm, 
               statistic = theta_proporcao, 
               R = 4000)

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

glimpse(ci)
## Observations: 1
## Variables: 5
## $ statistic <dbl> 0.242366
## $ bias      <dbl> 5.754127e-05
## $ std.error <dbl> 0.006079474
## $ conf.low  <dbl> 0.2301102
## $ conf.high <dbl> 0.2541969

Para a proporção temos o intervalo [-0.1844012, 0.06898568]

theta_correlacao = function(df,i) {
  df %>%
    slice(i) %>%
    summarise(correlacao = cor(mediana_pop[mediana_pop > 5], proporcao[mediana_pop > 5], method = "pearson" ) ) %>% pull(correlacao)
}

booted <- boot(data = lastfm, 
               statistic = theta_correlacao, 
               R = 4000)

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

glimpse(ci)
## Observations: 1
## Variables: 5
## $ statistic <dbl> -0.05679804
## $ bias      <dbl> 0.0005494606
## $ std.error <dbl> 0.06370502
## $ conf.low  <dbl> -0.1828482
## $ conf.high <dbl> 0.06422574