Implementando ICs - Lab 2.3 - FPCC2

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

Respondendo as Questões

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

Como solicitado, as métricas devem ser calculadas com base em uma amostras de tamanho igual a 300 de um universo populacional de 11.989 usuários disponível em nosso base de dados. É preciso primeiramente selecionar de maneira aleatoria essa amostra (n = 300). Em seguida, vamos construir o intervalo de confiança (IC) desta amostra para estimar a faixa de valores (mínimo e máximo) que a variável news apresenta em comparação a população.

Vamos encontrar a proporção de novos artistas em geral ter escutado por usuário é igual à quantidade de artivas novos dividido pela (quantidade artitas novos + quantidade artitas antigos). news/(news + old).

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

## Observations: 300
## Variables: 4
## $ 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...
## $ proportion  <dbl> 0.31460674, 0.15283843, 0.15662651, 0.20000000, 0.09722...

lastfm %>% ggplot(aes(x = proportion)) +
  geom_histogram(binwidth = .1,
                 colour = "darkred",
                 fill = "white")+ 
  labs(title = "Distribuição da proporção", x='Novos artistas', y='Usuários')

Como é possível verificar neste histograma, os valores da amostra apresentam uma proporção de novos artistas escutados próximo de 25% (0.25). Mas, não podemos tirar conclusões finais em apenas uma única amostra, é preciso fazer mais calculos. Vamos agora cálcular a média da proporção desta amostra:

Calculando pelo Bootstrapping manual

funcao_theta = function(df) {
  df %>%
    pull(proportion) %>%
    mean()
}
theta_c = funcao_theta(lastfm)

print(paste('Valor da Média amostral:',theta_c))

## [1] "Valor da Média amostral: 0.248356815094237"

Agora temos o valor de 0.248356815094237 para a média amostral. Vamos agora refazer o calculo para outras amostrar e calcular o Intervalo de Confiança. É necessário encontrar o número adequado de amostrar, o estudo indica um número de 5000 amostrar.

Calculando todas as novas amostras

set.seed(4321)

iteration = 5000

um_bootstrap <- function(x){
  proportion = x %>% pull(proportion)
  boot_x <- sample(proportion,           
                   size = NROW(proportion), 
                   replace = TRUE)
  return(mean(boot_x))
}
# A REAMOSTRAGEM
resampling = data.frame(i = 1:iteration) %>% 
  mutate(theta_c_s = map_dbl(i, ~ um_bootstrap(lastfm)))

Identificando o erro médio por amostragem

resampling_mean_error = resampling %>% 
  mutate(error = theta_c_s - theta_c)

resampling_mean_error

Calculando os Intervalos de Confiança

confidence_level = .95
alpha = 1 - confidence_level

interval = resampling_mean_error %>% 
  summarise(error_i = quantile(error, alpha/2), 
            error_s = quantile(error, 1 - alpha/2)) %>% 
  mutate(valor_i = theta_c + error_i, 
         valor_s = theta_c + error_s)

cat(paste('Est. iuperior:',interval$valor_i),"\n",
    paste('Est. snferior:',interval$valor_s),"\n",
    paste('Erro1:',interval$error_s),"\n",
    paste('Erro2:',interval$error_i)
    )

## Est. iuperior: 0.235798185658249 
##  Est. snferior: 0.261798064261403 
##  Erro1: 0.0134412491671663 
##  Erro2: -0.0125586294359884

Calculando pelo Bootstrapping com a biblioteca boot

set.seed(45678)
theta <- function(df,i) {
    mean(
        (df %>%
        slice(i) %>%
        pull(proportion)
    ))
}
booted <- boot(data = lastfm, 
               statistic = theta, 
               R = 5000)
ci = tidy(booted, 
          conf.level = .95,
          conf.method = "bca",
          conf.int = TRUE)

cat(paste('Média:',ci$statistic),"\n",
    paste('Est. inferior:',ci$conf.low),"\n",
    paste('Est. superior:',ci$conf.high),"\n",
    paste('Erro1:',ci$std.error), "\n",
    paste('Erro2:',ci$bias)
    )

## Média: 0.248356815094237 
##  Est. inferior: 0.235224000231867 
##  Est. superior: 0.262005093611465 
##  Erro1: 0.00676914569931786 
##  Erro2: -9.02014927335049e-05

Calculando os Intervalos de Confiança

ic_amostras = resampling_mean_error %>% 
  mutate(interval_i = theta_c_s + interval$error_i, 
         interval_s = theta_c_s + interval$error_s) %>% 
  mutate(contem_theta = theta_c >= interval_i & theta_c <= interval_s) 

ic_amostras %>% 
  summarise(cobertura = sum(contem_theta) / n())

ic_amostras %>% 
  slice(1:100) %>% 
  ggplot(aes(
    x = i,
    y = theta_c_s,
    ymin = interval_i,
    ymax = interval_s,
    color = contem_theta
  )) +
  geom_pointrange(alpha = .8, size = .3) +
  geom_hline(yintercept = theta_c, color = "dark blue") +
  labs(x = "Amostra",
       y = "Média") +
  scale_color_manual(values = c("red", "blue"))

Nesta representação gráfica dos valores dos IC que foram calculados podemos concluir que para 100 novas amostras, 5 delas não estão inseridas na média da população, atendendo ao nível de confiança de 95%.

Conclusão

É possível verificar que a média para ambas as bibliotecas foram iguais: 0.248356815094237. Agora quando verificamos os IC, vamos encontrar: o primeiro foi de [0.2357982,0.2617981] e o segundo foi de [0.235224,0.2620051]. Assim, podemos garantir que em média a proporção para novos artistas escutados pelos usuários é de 0.242366 com 95% de confiança e com um IC de [0.2357982,0.2617981].

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.

Precisamos, inicialmente, filtrar os usuários na amostra que possuem mediana_pop maior que 5. Após isso, podemos usar a correlação como a variável do intervalo de confiança. Iremos usar Spearman como o método de cálculo da correlação, por se tratar de uma amostra consideravelmente grande.Em seguida vamos utilizar a função bootstrap, que utiliza uma tabela para coletar as amostras de mesmo tamanho amostral.

Analisando a correlação entre as variáveis

lastfm %>%
    filter(mediana_pop >5) %>%
    mutate(ratio = news/(news+old)) %>%
    ggplot(aes(x=mediana_pop, y=ratio)) +
    geom_point()+
    geom_smooth(method='lm', color='red', se = FALSE)+
    labs(title='Correlação das variáveis', x = 'Gosto por música pop', y = 'Proporção novos artistas')

## `geom_smooth()` using formula 'y ~ x'

Diante do gráfico acima, podemos verificar que existe uma relação fraca e negativa entre as variáveis (mediana_pop, ratio). Vale destacar que pode existir uma relação não linear.

Calculando o bootstrapping manual

Vamos agora calcular a correlação amostrar pelo método Spearman.

funcao_theta_pop = function(df) {
  x <- df %>%
    filter(mediana_pop >5) %>%
    mutate(ratio = news/(news+old)) %>%
    summarise(corr = cor(mediana_pop, ratio, method= 'spearman')) %>%
    pull(corr) 
  return(x)
}
theta_c_pop = funcao_theta_pop(lastfm)

print(paste("Correlação amostral (spearman)",theta_c_pop))

## [1] "Correlação amostral (spearman) -0.0239004973681646"

Calculando as novas amostragem

repeticoes = 5000 

um_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_c_s_pop = map_dbl(i, ~ um_bootstrap(lastfm)))

Calculando o erro por amostragem

reamostragens_error = reamostragens %>% 
  mutate(error = theta_c_s_pop - theta_c_pop)

reamostragens_error

Calculando o intervalo de confiança

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

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

cat(paste('Est. iuperior:',interval$valor_i),"\n",
    paste('Est. snferior:',interval$valor_s),"\n",
    paste('Erro1:',interval$erro_i),"\n",
    paste('Erro2:',interval$erro_s)
    )

## Est. iuperior: 0.235798185658249 
##  Est. snferior: 0.261798064261403 
##  Erro1:  
##  Erro2:

É possível verificar que 95% dos valors do intervalo de confiança estão entre -0.146721481616017 e 0.102676021391378.

Calculando o bootstrapping pela biblioteca boot

theta_pop <- function(df,i) {
     df = df %>%
        filter(mediana_pop > 5) %>%
        slice(i) %>%
        mutate(prop = news/(news + old),
               cor = cor(mediana_pop, prop, method= 'spearman'))
    mean(df$cor)
}
booted_pop <- boot(data = lastfm, 
               statistic = theta_pop, 
               R = 5000)
ci = tidy(booted_pop, 
          conf.level = .95,
          conf.method = "bca",
          conf.int = TRUE)

cat(paste('Correlação amostral:',ci$statistic),"\n",
    paste('Est. inferior:',ci$conf.low),"\n",
    paste('Est. Superior:',ci$conf.high),"\n",
    paste('Erro1:',ci$std.error), "\n",
    paste('Erro2:',ci$bias)
    )

## Correlação amostral: -0.0239004973681646 
##  Est. inferior: -0.151558294004748 
##  Est. Superior: 0.100900483580033 
##  Erro1: 0.0639430535889186 
##  Erro2: 0.000501077366240138

Agora com a utilização da biblioteca boot não vamos encontrar uma igualdade entre os resultados, lembrando que a cada execução de novos valores da amostrar os resultados são diferentes e quanto maior o número de novas amostras mais resultados tendem a se aproximar.

Conclusão

É possível verificar que o valor do theta do bootstrapping manual foi de -0.0239004973681646 e pela biblioteca boot foi -0.0239004973681646, ambos iguais, já o IC foram de [-0.146721481616017, 0.102676021391378] e no segundo [-0.149875704194716,0.0991188594145701], respectivamente, agora com valores bem próximos. Com 95% de confiança é possível confirmar que, utilizando o método de spearman, a correlação é de -0.0239004973681646 entre a popularidade mediana dos artistas escutado e a proporção dos artistas escutados que eram novos.