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

Crie intervalos com 95% de confiança.

Proporção geral de artistas novos escutados

A proporção de novos artistas entre os ouvidos será calculada através da razão entre a quantidade de artistas novos e o total de artistas ouvidos por um usuário (news / news + old).

proporcoes <- lastfm %>% 
  mutate(proporcao = news / (news + old))
proporcoes %>% 
  ggplot(aes(proporcao)) +
  geom_histogram(binwidth = 0.05,
                 fill = "white",
                 colour = "darkgrey") +
  labs(title = "Ditribuição das proporções de novos artistas ouvidos")

Dos artistas mais ouvidos pela maioria dos usuários analisados, 20% são novos. Utilizaremos a média para generalizar a proporção de novos artistas ouvidas pelos usuários do LastFM.

Bootstraping manual

funcao_theta_prop <- function(df) {
  media <- df %>% 
    pull(proporcao) %>% 
    mean()
  return(media)
}

theta_prop <- funcao_theta_prop(proporcoes)
repeticoes = 4000

prop_bootstrap <- function(data) {
  props <- data %>% pull(proporcao)
  boot_data <- sample(props,
                      size = NROW(proporcoes),
                      replace = TRUE)
  media <- mean(boot_data)
  return(media)
}

set.seed(1212)

reamostragens <- tibble(i = 1:repeticoes) %>% 
  mutate(theta_prop_boot_m = map_dbl(i, ~ prop_bootstrap(proporcoes)))

reamostragens <- reamostragens %>% 
  mutate(erro = theta_prop_boot_m - theta_prop)

Vejamos a distribuição amostral da proporção e a do erro entre a proporção da nossa amostra original e as reamostragens do bootstraping. Nos gráficos abaixo, a linha preta na vertical indica o valor da proporção média para a amostra original e o valor zero, respectivamente. Utilizaremos o zero como valor ideal para a diferença, pois, nestas situações, a estatística para uma reamostragem foi igual à estatística da amostra original.

reamostragens %>% 
  ggplot(aes(theta_prop_boot_m)) +
  geom_histogram(binwidth = .0025,
                 fill = "white",
                 colour = "darkgrey") +
  geom_vline(xintercept = theta_prop) +
  labs(title = "Distribuição amostral",
       x = "Proporção média")

Vemos que o valor da proporção média não difere muito da média original em boa parte das reamostragens.

reamostragens %>% 
  ggplot(aes(erro)) +
  geom_histogram(binwidth = .0025,
                 fill = "white",
                 colour = "darkgrey") +
  geom_vline(xintercept = 0) +
  labs(title = "Distribuição do erro amostral")

Analisando os erros, podemos observar melhor que o valor geral da proporção para a maioria das amostras não difere muito do da amostra original.

Agora podemos calcular o intervalo de confiança para a média da proporção de novos artistas ouvidos.

intervalo = reamostragens %>% 
  summarise(erro_i = quantile(erro, .05), 
            erro_s = quantile(erro, .95)) %>% 
  mutate(valor_i = theta_prop + erro_i, 
         valor_s = theta_prop + erro_s)
intervalo
ggplot() +
  geom_rect(data = intervalo,
            aes(xmin = valor_i, xmax = valor_s),
            ymin = -Inf,
            ymax = Inf,
            fill = "red",
            alpha = .25) +
  geom_histogram(data = reamostragens,
                 aes(x = theta_prop_boot_m),
                 binwidth = 0.0025,
                 fill = "white",
                 colour = "darkgrey") +
  geom_vline(xintercept = theta_prop) +
  labs(title = "Intervalo ao redor da média original",
       x = "Proporção média")

Bootstraping com a biblioteca boot

Utilizaremos a biblioteca boot para calcular o intervalo de confiança para a proporção média de novos artistas ouvidos.

repeticoes = 4000

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

booted <- boot(data = proporcoes, 
               statistic = funcao_theta_prop_lib, 
               R = repeticoes)

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

ci

O valor da proporção geral de novos artistas ouvidos pelos usuários foi de 0.242366 em ambas as abordagens. O intervalo de confiança na abordagem manual foi [0.2325967, 0.2524944] e utilizando a biblioteca boot foi de [0.2311165, 0.2546783]. Com 95% de confiança, podemos dizer que a proporção média de novos artistas ouvidos pelos usuários do LastFM é de 0.242366, ou seja, aproximadamente 20% dos artistas ouvidos são novos.

ggplot() +
  geom_errorbar(aes(x = 'Manual',
                    ymin = intervalo$valor_i,
                    ymax = intervalo$valor_s),
                width = 0.095) +
  geom_point(aes(x = 'Manual',
                 y = theta_prop)) +
  geom_errorbar(aes(x = 'Biblioteca Boot',
                    ymin = ci$conf.low,
                    ymax = ci$conf.high),
                width = 0.095) +
  geom_point(aes(x = 'Biblioteca Boot',
                 y = theta_prop)) +
  labs(title = "Intervalos gerados em cada abordagem",
       x = NULL,
       y = "Proporção média")

Relação entre popularidade e proporção

Vejamos inicialmente como ser comporta a proporção em função da popularidade.

mais_pops <- proporcoes %>% 
  filter(mediana_pop > 5)
mais_pops %>% 
  ggplot(aes(x = mediana_pop,
             y = proporcao)) +
  geom_point() +
  labs(title = "Proporção em função da popularidade",
       x = "Popularidade",
       y = "Propoção")

Não é possível identificar uma relação clara entre as duas variáveis. Usaremos o valor de correlação como estatística para quantificar a relação entre a proporção e popularidade.

Bootstraping manual

funcao_theta_pop <- function(df){
  x <- df %>% 
    summarise(corr = cor(mediana_pop, proporcao, method = 'pearson')) %>% 
    pull(corr)
  return(x)
}

theta_corr <- funcao_theta_pop(mais_pops)
repeticoes = 4000

corr_bootstrap <- function(df){
  boot_x <- sample_n(df,           
                   size = NROW(df),
                   replace = TRUE) 
  return(funcao_theta_pop(boot_x))
}

reamostragens <- tibble(i = 1:repeticoes) %>% 
  mutate(theta_corr_boot_m = map_dbl(i, ~ corr_bootstrap(mais_pops)))

reamostragens <- reamostragens %>% 
  mutate(erro = theta_corr_boot_m - theta_corr)
reamostragens %>% 
  ggplot(aes(theta_corr_boot_m)) +
  geom_histogram(binwidth = .025,
                 fill = "white",
                 colour = "darkgrey") +
  geom_vline(xintercept = theta_corr) +
  labs(title = "Distribuição amostral",
       x = "Correlação")

Pelo gráfico vemos que para a maioria das amostras o valor da correlação é tão baixo quanto para a amostra original.

reamostragens %>% 
  ggplot(aes(erro)) +
  geom_histogram(binwidth = .025,
                 fill = "white",
                 colour = "darkgrey") +
  geom_vline(xintercept = 0) +
  labs(title = "Distribuição do erro amostral")

Analisando os erros vemos que a maioria das amostras obteve o mesmo valor de correlação que o calculado na amostra original.

intervalo <- reamostragens %>% 
  summarise(erro_i = quantile(erro, .05), 
            erro_s = quantile(erro, .95)) %>% 
  mutate(valor_i = theta_corr + erro_i, 
         valor_s = theta_corr + erro_s)
ggplot() +
  geom_rect(data = intervalo,
            aes(xmin = intervalo$valor_i, 
                xmax = intervalo$valor_s),
            ymin = -Inf,
            ymax = Inf,
            fill = "red",
            alpha = .25) +
  geom_histogram(data = reamostragens,
                 aes(x = theta_corr_boot_m),
                 binwidth = 0.025,
                 fill = "white",
                 colour = "darkgrey") +
  geom_vline(xintercept = theta_corr) +
  labs(title = "Intervalo ao redor da correlação original")

Bootstraping com a biblioteca boot

Agora utilizando a biblioteca boot.

repeticoes = 4000 

funcao_theta_corr_lib <- function(df, i){
  x <- df %>% 
    slice(i) %>% 
    summarise(correlacao = cor(mediana_pop, proporcao, method = 'pearson')) %>% 
    pull(correlacao)
  return(x)
}

booted <- boot(data = mais_pops, 
               statistic = funcao_theta_corr_lib, 
               R = repeticoes)

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

ci

A correlação entre proporção e popularidade calculada foi de -0.056798, o que indica uma relação muito fraca entre as variáveis. O intervalo de confiança calculado na abordagem manual foi de [-0.157639, intervalo %>% pull(valor_s)] e utilizando a biblioteca boot foi de [-0.1841781, 0.064871]. Com 95% de confiança, podemos dizer que a relação entre a proporção de artistas novos ouvidos e a mediana da popularidade dos artistas ouvidos por um usuário é muito baixa.

ggplot() +
  geom_errorbar(aes(x = 'Manual',
                    ymin = intervalo$valor_i,
                    ymax = intervalo$valor_s),
                width = 0.095) +
  geom_point(aes(x = 'Manual',
                 y = theta_corr)) +
  geom_errorbar(aes(x = 'Biblioteca Boot',
                    ymin = ci$conf.low,
                    ymax = ci$conf.high),
                width = 0.095) +
  geom_point(aes(x = 'Biblioteca Boot',
                 y = theta_corr)) +
  labs(title = "Intervalos gerados em cada abordagem",
       x = NULL,
       y = "Correlação")