Os dados

lastfm = read_csv(here("dados/experimento-lastfm.csv"), 
                  col_types = cols(.default = col_double(), 
                                   user = col_character()))
glimpse(lastfm)
## Observations: 11,989
## Variables: 15
## $ user           <chr> "00blixt", "0ctane", "0sebastian0", "0xkirill", "

## $ ecletic        <dbl> 1159.4960, 3439.8373, 1611.4464, 697.0974, 1660.1

## $ clusters_total <dbl> 37, 155, 68, 37, 72, 93, 163, 150, 88, 132, 130, 

## $ clusters       <dbl> 19, 57, 25, 13, 31, 37, 64, 59, 27, 58, 62, 12, 3

## $ media_pop      <dbl> 5.536362, 4.726167, 5.804086, 4.997963, 5.800357,

## $ mediana_pop    <dbl> 5.660406, 4.730935, 5.896073, 5.214290, 5.865121,

## $ dp_pop         <dbl> 0.6852190, 0.8819529, 0.6167104, 0.8685811, 0.544

## $ news           <dbl> 38, 32, 24, 36, 28, 25, 30, 28, 14, 22, 27, 21, 2

## $ news_total     <dbl> 38, 32, 24, 36, 28, 25, 30, 28, 14, 22, 27, 21, 2

## $ old            <dbl> 116, 130, 10, 46, 128, 95, 228, 53, 46, 118, 90, 

## $ old_total      <dbl> 119, 141, 13, 71, 133, 167, 256, 63, 52, 125, 102

## $ count_sim      <dbl> 0.155514371, 0.074823398, -0.123188406, 0.0841938

## $ tempo_sim      <dbl> 0.208395596, -0.020325875, -0.043478261, -0.00160

## $ count_pop      <dbl> 0.26923104, 0.13144651, -0.16666667, 0.09372521, 

## $ tempo_pop      <dbl> 0.08303680, -0.14228113, -0.11594203, 0.13012993,

lastfm = select(lastfm, news, ecletic)

lastfm %>% ggplot(aes(news)) + geom_histogram(binwidth = 10)

lastfm %>% ggplot(aes(ecletic)) + geom_histogram(binwidth = 100)

Amostra e visÔes a partir dela

Imagine por agora que os dados que temos do Last.fm são completos. Que são a população inteira de usuårios. Como seria nossa visão dos dados se tivéssemos apenas uma amostra dos dados? sample(x, n) faz uma amostra aleatórioa de n elementos tirados do vetor x:

Se calcularmos a mĂ©dia do nĂșmeros de novos artistas escutados para trĂȘs amostras de 100 elementos, teremos 3 resultados diferentes:

lastfm %>% pull(news) %>% sample(100) %>% mean()
## [1] 31.95
lastfm %>% pull(news) %>% sample(100) %>% mean()
## [1] 35.37
lastfm %>% pull(news) %>% sample(100) %>% mean()
## [1] 34.02

Se fizermos isso muitas vezes podemos ver como essa variação acontece. A distribuição dos valores de uma estatística em diferentes amostras de uma população se chama distribuição amostral da estatística.

set.seed(1)

amostras = data.frame(amostra = 1:1000) %>% # faremos 1000 vezes
  mutate(media = map_dbl(amostra, ~ lastfm %>% 
                           pull(news) %>% 
                           sample(100) %>% 
                           mean()))

amostras
amostras %>%
  ggplot(aes(media)) +
  geom_histogram(binwidth = .5,
                 fill = "white",
                 colour = "darkgrey") +
  geom_vline(xintercept = mean(lastfm$news)) # theta: média calculada com todos os dados

E se o tamanho da amostra (n) fosse muito menor?

amostras = data.frame(amostra = 1:1000) %>% # faremos 1000 vezes
  mutate(media = map_dbl(amostra, ~ lastfm %>% 
                                      pull(news) %>% 
                                      sample(10) %>% 
                                      mean()))

amostras
amostras %>% 
  ggplot(aes(media)) + 
  geom_histogram(binwidth = .5, fill = "white", colour = "darkgrey") + 
  geom_vline(xintercept = mean(lastfm$news)) # média calculada com todos os dados

A mesma lógica vale para outras estatísticas além da média. O código abaixo analisa a distribuição dos valores observados em amostras a partir das quais calculamos a mediana. Altere o código para usar outra estatística: (dica: max e min não funcionam.)

amostras = data.frame(amostra = 1:1000) %>% # faremos 1000 vezes
  mutate(estatistica = map_dbl(amostra, ~ lastfm %>% 
                                      pull(news) %>% 
                                      sample(size = 100) %>% 
                                      median()))

amostras
amostras %>% 
  ggplot(aes(estatistica)) + 
  geom_histogram(binwidth = 1, fill = "white", colour = "darkgrey") + 
  geom_vline(xintercept = median(lastfm$news)) # média calculada com todos os dados

Distribuição amostral

A distribuição dos valores de uma estatística a partir de amostras é chamada de distribuição amostral daquela estatística. Ela tem um papel importante, porque é a partir do entendimento dela que estimaremos quanta confiança temos em uma estatística que estamos calculando a partir de uma amostra.

O principal a entender aqui é que se conhecermos a distribuição amostral, saberemos quão longe normalmente a estatística calculada para uma amostra estå daquela calculada para a população. Sabendo isso, podemos calcular uma margem de erro para a estimativa feita a partir da amostra.

A situação mais comum na vida real

Normalmente temos uma amostra da população apenas. Daí não conhecemos a distribuição amostral. Mas gostaríamos de a partir da nossa amostra estimar onde estå a estatística para a população.

Exemplo: queremos estimar qual a proporção de pessoas que gostarå do produto (a estatística) para o universo de todos os usuårios (a população) usando o cålculo da proporção de pessoas que gostou do produto (a mesma estatística) em um teste com 100 pessoas (a amostra).

Se conhecermos como a estatística varia na distribuição amostral (ex: 2 pontos pra mais ou pra menos cobrem 99% dos casos) e temos a estatística calculada para a amostra, poderíamos estimar uma faixa de valores onde achamos que a estatística estå para a população com 99% de confiança.

A ideia central que usaremos

Não temos a população, para estimar a variação na distribuição amostral. A ideia principal que usaremos, em uma técnica chamada boostrapping é que usar a amostra como substituto da população e simular a amostragem através de reamostragem com reposição fornece uma estimativa precisa da variação na distribuição amostral.

Pegando por partes:

  • Consideramos a amostra \(A\) que tem tamanho \(n\) como sendo um substituto da população
  • Repetimos \(b\) vezes o seguinte processo: criamos uma amostra de tamanho \(n\) obtendo elementos aleatoriamente de \(A\), repondo cada elemento depois de cada sorteio.
  • Calculamos a estatĂ­stica \(e\) que nos interessa (mĂ©dia, mediana, desvio padrĂŁo, o que for) para cada uma das \(b\) amostras.

Como resultado, sabemos como a estatística \(e\) varia em uma simulação de \(b\) amostragens.

O princĂ­pio do bootstrap diz que a variação de \(e\) nos bootstraps aproxima a variação de \(e_a\) em mĂșltiplas amostras.

Aplicando bootstrapping

Conhecendo a variação de \(e_a\) em torno de \(e_p\) podemos estimar quão longe a estatística \(e_a\) que calculamos para nossa amostra provavelmente estå da estatística \(e_p\) para a população. Isso porque se sabemos normalmente quão longe \(e_a\) estå de \(e_p\) (que é o que estimamos com o boostrap), sabemos também o inverso.

Por exemplo, se estimamos que em 99% das amostras, \(e_a\) estå a +/- 10 unidades de \(e_p\), e calculamos que \(e_a = 22\), podemos estimar que _com 99% de confiança, \(e_p\) estå em \(e_a +/- 10 = [12, 32]\).

No nosso caso nós não temos a variação de \(e_a\). Porém contornamos essa dificuldade via o princípio do bootstrapping, que diz que a distribuição de \(e_b\) em muitos bootstraps aproxima a variação de \(e_a\).

Caso conhecĂȘssemos a distribuição de \(e_a\), e definĂ­ssemos \(\delta = e_a - e_p\), para calcular um intervalo onde essa diferença cai \(.90\) do tempo nĂłs conseguirĂ­amos achar os percentis \(\delta_{0.05}\) e \(\delta_{0.95}\) tal que:

\(P(\delta_{0.05} \le e_a - e_p \le \delta_{0.95} | e_p) = .9\). Manipulando, temos que \(P(e_a - \delta_{0.05} \ge e_p \ge e_a - \delta_{0.95} | e_p) = .9\). Ou seja, o IC Ă© \([e_a - \delta_{0.05}; e_a - \delta_{0.95}]\).

Pelo princípio do boostrap, podemos aplicar a mesma lógica, porém usando a distribuição do bootstrap para substituir a amostral.

Aplicando bootstrapping

amostra = lastfm %>%  
  sample_n(200) %>% 
  pull(news)

media_amostra = mean(amostra)

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

um_bootstrap <- function(x){
  boot_x <- sample(x,
                   size = NROW(x), 
                   replace = TRUE) # aqui Ă© o bootstrap
  return(mean(boot_x))
}

# Agora o input das reamostragens Ă© a amostra!
experimentos = tibble(i = 1:repeticoes) %>% 
  group_by(i) %>% 
  mutate(media_bootstrap = map_dbl(i, ~ um_bootstrap(amostra)))

ggplot(experimentos, aes(x = media_bootstrap)) + 
  geom_histogram(binwidth = 1, colour = "darkorange", fill = "white")

summary(experimentos)
##        i        media_bootstrap
##  Min.   :   1   Min.   :25.56  
##  1st Qu.:1001   1st Qu.:30.73  
##  Median :2000   Median :32.05  
##  Mean   :2000   Mean   :32.10  
##  3rd Qu.:3000   3rd Qu.:33.35  
##  Max.   :4000   Max.   :39.38

Calculando o IC:

# IC com 90%: 
alpha = .1
quantile(experimentos$media_bootstrap, 
         probs = c(.05, .95))
##       5%      95% 
## 29.11500 35.39025
cis = experimentos %>% 
  ungroup() %>% 
  mutate(diferenca = media_bootstrap - media_amostra) %>% 
  summarise(l = quantile(diferenca, alpha/2), 
            u = quantile(diferenca, 1 - alpha/2)) %>% 
  mutate(ci_lower = media_amostra - u, 
         ci_upper = media_amostra - l)

ggplot(experimentos, aes(x = media_bootstrap)) + 
  geom_histogram(binwidth = 1, colour = "darkorange", fill = "white") + 
  geom_vline(aes(xintercept = cis$ci_lower), colour = "blue", size = 1.5) + 
  geom_vline(aes(xintercept = cis$ci_upper), colour = "blue", size = 1.5) + 
  geom_vline(aes(xintercept = media_amostra), colour = "green") + 
  geom_vline(aes(xintercept = mean(lastfm$news)), colour = "red")

Com outro nível de confiança:

# IC com 95%: 
alpha2 = .05
cis2 = experimentos %>% 
    mutate(diferenca = media_bootstrap - media_amostra) %>% 
    summarise(l = quantile(diferenca, alpha2/2), 
              u = quantile(diferenca, 1 - alpha2/2)) %>% 
    mutate(ci_lower = media_amostra - u, 
           ci_upper = media_amostra - l)

ggplot(experimentos, aes(x = media_bootstrap)) + 
  geom_histogram(binwidth = 1, colour = "pink", fill = "white") + 
    geom_vline(aes(xintercept = cis$ci_lower), colour = "blue") + 
  geom_vline(aes(xintercept = cis$ci_upper), colour = "blue") + 
  geom_vline(aes(xintercept = cis2$ci_lower), colour = "purple", size = 2) + 
  geom_vline(aes(xintercept = cis2$ci_upper), colour = "purple", size = 2) + 
  geom_vline(aes(xintercept = media_amostra), colour = "green")