Implementando ICs

Os dados

set.seed(12345)

read_csv(here::here("data/experimento-lastfm.csv"),
         col_types = cols(.default = col_double(),
                          user = col_character()))%>%
  sample_n(300) %>% 
  select(news, old, mediana_pop) %>%
  mutate(prop_news =  news / (news + old)) -> lastfm

lastfm %>%
  glimpse()

## Observations: 300
## Variables: 4
## $ 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…
## $ prop_news   <dbl> 0.19531250, 0.17567568, 0.25301205, 0.09090909, 0.1538461…

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.

Proporção de novos artistas em geral

Observermos nossa amostra logo abaixo através de um gráfico de dispersão

lastfm %>%
  ggplot(aes("", prop_news)) +
  geom_jitter() +
  stat_summary(fun = median,
               fun.min = median, fun.max = median,
               geom = "crossbar", color="darkred",
               width = 0.5) +
  labs(x= "", y="Proporção de artistas novos")

• Mediana próxima de 0.2. Isso aponta que para metade dos usuários a proporção de artistas novos que estes escutam chega até 0.2 (20% da mídia consumida é de artistas novos).
• Vemos um usuário em particular para o qual a proporção de artistas novos foi de 0.6 (60% dos artistas foi de artistas novos).

Aplicando biblioteca boot

my_theta <- function(x, i) {
  x %>% 
    slice(i) %>% 
    summarise(median_prop_news = median(prop_news)) %$%
    median_prop_news -> result
    
    return(result)
}

median_prop.boot <- boot(data = lastfm, 
                 statistic = my_theta, 
                 R = 5000)

plot(median_prop.boot)

• Avaliação dos quantis e da distribuição da estatística gerada via reamostragem não sugerem forte viés ou assimetria que desencorajariam o uso do “basic” bootstrap.

boot.ci(boot.out = median_prop.boot, 
        conf = 0.95, 
        type = "basic") -> median_prop.ci

median_prop.ci

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 5000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = median_prop.boot, conf = 0.95, type = "basic")
## 
## Intervals : 
## Level      Basic         
## 95%   ( 0.2143,  0.2452 )  
## Calculations and Intervals on Original Scale

median_prop.data <- tidy(median_prop.boot)
median_prop.data$conf.low = median_prop.ci$basic[4]
median_prop.data$conf.high = median_prop.ci$basic[5]
median_prop.data$origem <- "Biblioteca boot"

median_prop.data

• A mediana da amostra inicial de fato está próxima de 0.2 (0.2280154)
• Baixa variabilidade e bias (Pouca incerteza a ser representada no intervalo de confiança)
• O intervalo de confiança é de tamanho relativamente pequeno [0.2142727, 0.2451875
  ].

Através do intervalo de confiança temos uma representação próxima da apontada pela amostra inicial, estimamos que em geral pouco mais de 20% (proporção próxima de 0.2) do que os usuários consomeme se refere a artistas novos

Implementação própria

Ao longo deste relatório utilizaremos uma implementação própria do “basic bootstrap”

lastfm %>%
  summarise(median_prop_news = median(prop_news)) %$%
  median_prop_news -> r


estimates <- c()

estimates <- replicate(10000, median(sample(lastfm$prop_news, replace=TRUE)))

q <- quantile(estimates, probs=c(0.025,0.975))
lower <- 2 * r - q[2]
upper <- 2 * r - q[1]

cat("Intervalo de confiança: \n", "[",lower,",",upper,"]", sep="")

## Intervalo de confiança: 
## [0.2141614,0.2447406]

• Existe uma divergência pequena entre os resultados de nossa implementação e os resultados da biblioteca. A diferença é mais alta no limite superior

Apliquemos nosso método de cálculo à reamostragem provida pela biblioteca boot :

q <- quantile(median_prop.boot$t, probs=c(0.025,0.975))

lower <- 2 * r - q[2]
upper <- 2 * r - q[1]

cat("Intervalo de confiança: \n", "[",lower,",",upper,"]", sep="")

## Intervalo de confiança: 
## [0.2142727,0.2451875]

• Os resultados batem com os da biblioteca, estamos de fato calculando o intervalo de confiança corretamente. As divergências partem de particularidades na maneira com que a biblioteca boot executa sua reamostragem.

Coloquemos nossos resultados em formato de dataframe pra facilitar a criação de visualizações

conf.low = c(lower)
conf.high = c(upper)
inhouse_median_prop.data <- data.frame(conf.low,conf.high)
inhouse_median_prop.data$statistic <- r
inhouse_median_prop.data$bias <- NA
inhouse_median_prop.data$std.error <- NA

inhouse_median_prop.data$origem <- "Nossa implementação"

inhouse_median_prop.data

rbind(inhouse_median_prop.data,
      median_prop.data) -> median_data

median_data

median_data %>%
  ggplot(aes(x = origem, y=statistic,
             ymin = conf.low, ymax = conf.high)) +
  geom_errorbar(width = .15) +
  geom_point(color="red", size=2) +
  scale_y_continuous(limits = c(0, 1)) +
  labs(x="", y="Proporção de artistas novos",
       title=expression(paste("basic bootstrap ", alpha, " = 0.05")))

* Os resultados da biblioteca boot e de nossa implementação estão muito próximos.

Através do intervalo de confiança temos uma representação próxima da apontada pela amostra inicial, estimamos que em geral pouco mais de 20% (proporção próxima de 0.2) do que os usuários consomem se refere a artistas novos

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.

Significado do efeito

• Uma correlação positiva nos diria que existe uma tendência entre os usuários que escutam artistas mais populares a buscar mais novo conteúdo, ou uma tendência entre usuários que mais buscam conteúdo novo a escutar artistas mais populares.

• Uma correlação negativa nos diria que existe uma tendência entre os usuários que escutam artistas mais populares a se abster de buscar conteúdo novo, ou uma tendência entre usuários que se abstem de buscar conteúdo novo a escutar artistas mais populares.

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

lastfm_pop %>% 
  glimpse()

## Observations: 270
## Variables: 4
## $ 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…
## $ prop_news   <dbl> 0.19531250, 0.17567568, 0.25301205, 0.09090909, 0.1538461…

Observermos nossa amostra logo abaixo através de um gráfico de dispersão

lastfm_pop %>%
  ggplot(aes(mediana_pop, prop_news)) +
  geom_jitter() +
  geom_smooth(method='loess', formula='y ~ x') +
  labs(x= "Popularidade mediana dos artistas escutados",
       y="Proporção de artistas novos")

• De antemão não existe uma tendência clara de correlação na amostra inicial

Aplicando biblioteca boot

my_theta <- function(x, i) {
  x %>% 
    slice(i) %>%
    summarise(cor_prop_pop = cor(prop_news, mediana_pop,
                               method="kendall")) %$%
    cor_prop_pop -> result
    
    return(result)
}


lastfm_pop %>%
  boot(statistic = my_theta,
       R = 10000) -> cor.boot 

plot(cor.boot)

• Avaliação dos quantis e da distribuição da estatística gerada via reamostragem não sugerem forte viés ou assimetria que desencorajariam o uso do “basic” bootstrap.

boot.ci(boot.out = cor.boot, 
        conf = 0.95, 
        type = "basic") -> cor.ci

cor.ci

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = cor.boot, conf = 0.95, type = "basic")
## 
## Intervals : 
## Level      Basic         
## 95%   (-0.0731,  0.0789 )  
## Calculations and Intervals on Original Scale

O intervalo de confiança produzido inclui o valor 0 e está restrito a valores muito baixos de correlação. Não foi possível estabelecer um efeito estatisticamente significativo. A partir desta amostra estimamos que pode não haver um efeito correlação entre a popularidade mediana dos artistas escutado e a proporção dos artistas escutados que eram novos, ou se ele existir, ele provavelmente é pequeno em qualquer das direções.

cor.data <- tidy(cor.boot)
cor.data$conf.low = cor.ci$basic[4]
cor.data$conf.high = cor.ci$basic[5]
cor.data$origem = "Biblioteca boot"

cor.data

• O valor da estatística correlação na amostra inicial (statistic) é também muito próximo de zero.

Aplicando implementação própria

lastfm_pop %>%
  summarise(cor_prop_pop = cor(prop_news, mediana_pop,
                               method="kendall")) %$%
    cor_prop_pop -> r


estimates <- c()

for (i in 1:1000) {

  # Draw a bootstrap sample.
  sample_data <- lastfm_pop[sample(nrow(lastfm_pop), replace=T),]

  # Estimate the statistic.
  sample_data %>%
    summarise(cor_prop_pop = cor(prop_news, mediana_pop,
                                 method="kendall")) %$%
      cor_prop_pop -> statistic

  # Save the bootstrap estimate of the statistic 
  estimates <- append(estimates, statistic)

}

q <- quantile(estimates, probs=c(0.025,0.975))
lower <- 2 * r - q[2]
upper <- 2 * r - q[1]

cat("[",round(lower, digits = 4),",",round(upper, digits = 4),"]", sep="")

## [-0.0718,0.0781]

• Existe uma pequena divergência entre os resultados de nossa implementação e os resultados da biblioteca. Embora a dvergência seja relativamente baixa comparada ao estudo anterior vemos um aumento. A diferença é maior no limite inferior.

Apliquemos nosso método de cálculo à reamostragem provida pela biblioteca boot :

q <- quantile(cor.boot$t, probs=c(0.025,0.975))

lower <- 2 * r - q[2]
upper <- 2 * r - q[1]

cat("[",round(lower, digits = 4),",",round(upper, digits = 4),"]", sep="")

## [-0.073,0.0789]

• Os resultados batem com os da biblioteca, estamos de fato calculando o intervalo de confiança corretamente. As divergências partem de particularidades na maneira com que a biblioteca boot executa sua reamostragem.

Coloquemos nossos resultados em formato de dataframe pra facilitar a criação de visualizações

conf.low = c(lower)
conf.high = c(upper)
inhouse_cor.data <- data.frame(conf.low,conf.high)
inhouse_cor.data$statistic <- r
inhouse_cor.data$bias <- NA
inhouse_cor.data$std.error <- NA

inhouse_cor.data$origem <- "Nossa implementação"

inhouse_cor.data

rbind(inhouse_cor.data,
      cor.data) -> correlation_data

correlation_data

correlation_data %>%
  ggplot(aes(x = origem, y=statistic,
             ymin = conf.low, ymax = conf.high)) +
  geom_errorbar(width = .15) +
  geom_hline(yintercept = 0, colour = "darkorange") +
  geom_point(color="red", size=2) +
  scale_y_continuous(limits = c(-1, 1)) +

  labs(x="", y="Correlação (Popularidade / Proporção novos)",   
       title=expression(paste("basic bootstrap ", alpha, " = 0.05")))

• Os resultados da biblioteca boot e de nossa implementação são similares

O intervalo de confiança produzido inclui o valor 0 e está restrito a valores muito baixos de correlação (muito próximos de 0). Não foi possível estabelecer um efeito estatisticamente significativo. A partir desta amostra estimamos que pode não haver um efeito correlação entre a popularidade mediana dos artistas escutado e a proporção dos artistas escutados que eram novos, ou se ele existir, ele provavelmente é pequeno em qualquer das direções.