Os dados

Usamos os dados coletados por Andryw Marques para o mestrado dele. Cada medição reflete os hábitos de uma pessoa usando o last.fm durante um semestre. Usaremos apenas um atributo que crio abaixo: o número (médio) de artistas escutados por mês pela pessoa.

lastfm = read_csv(
  here::here("data/experimento-lastfm.csv"),
  col_types = cols(.default = col_double(),
                   user = col_character())
)

lastfm = select(lastfm, news, old, ecletic) %>% filter(!is.na(news + old)) %>% mutate(artists = (news + old)/6)

glimpse(lastfm)
## Rows: 11,988
## Columns: 4
## $ news    <dbl> 38, 32, 24, 36, 28, 25, 30, 28, 14, 22, 27, 21, 23, 37, 25, 3…
## $ old     <dbl> 116, 130, 10, 46, 128, 95, 228, 53, 46, 118, 90, 34, 100, 109…
## $ ecletic <dbl> 1159.4960, 3439.8373, 1611.4464, 697.0974, 1660.1563, 1907.90…
## $ artists <dbl> 25.666667, 27.000000, 5.666667, 13.666667, 26.000000, 20.0000…
lastfm %>% ggplot(aes(artists)) + geom_histogram(binwidth = 5, boundary = 0)

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. Se o que nos interessa por exemplo é a média dessa população, ela é exatamente:

THETA = lastfm %>%  
  pull(artists) %>% 
  mean() # theta: média calculada com todos os dados

THETA
## [1] 22.14312

Como seria nossa visão dos dados se tivéssemos apenas uma amostra dos dados? sample_n(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 (3 \(\hat{\theta}\) diferentes):

lastfm %>% sample_n(size = 100) %>% pull(artists) %>%  mean()
## [1] 21.065
lastfm %>% sample_n(100) %>% pull(artists) %>% mean()
## [1] 22.005
lastfm %>% sample_n(100) %>% pull(artists) %>% mean()
## [1] 21.91333

Se fizermos isso muitas vezes vemos como essa variação de \(\hat{\theta}\) 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)

amostra_theta_c = function(df, n = 100) {
  df %>%
    sample_n(n) %>% 
    pull(artists) %>% 
    mean()
}

amostras = tibble(amostra = 1:1000) %>% # faremos 1000 vezes
  mutate(theta_c = map_dbl(amostra, ~ amostra_theta_c(lastfm)))

amostras
amostras %>%
  ggplot(aes(theta_c)) +
  geom_histogram(binwidth = .5,
                 fill = "white",
                 colour = "darkgrey") +
  geom_vline(xintercept = THETA, linetype = "dashed") + 
  labs(title = "Distribuição amostral")

Podemos ver também qual a distribuição do erro amostral:

amostras = amostras %>% 
  mutate(erro = theta_c - THETA)

amostras
amostras %>%
  ggplot(aes(erro)) +
  geom_histogram(binwidth = .5,
                 fill = "white",
                 colour = "darkblue") +
  geom_vline(xintercept = 0, linetype = "dashed") + 
  labs(title = "Distribuição do erro amostral")

Veja que conhecendo a distribuição de \(\delta = \hat{\theta} - \theta\) nós conhecemos como as visões de \(\theta\) a partir de amostras variam. Com isso nós sabemos quão longe \(\hat{\theta}\) geralmente está de \(\theta\) e vice-versa.

Na prática: usando a distribuição nós conseguimos encontrar 2 valores de \(\delta = \hat{\theta} - \theta\) entre os quais os quais \(\hat{\theta} - \theta\) está 90% do tempo. Basta encontrar o 5 e o 95 percentis. Esses são os dois valores mais próximos que contém 90% das observações.

intervalo = amostras %>% 
  summarise(erro_i = quantile(erro, .05), 
            erro_s = quantile(erro, .95))

intervalo

Esse esse intervalo ao redor de \(\theta\) é o intervalo onde \(\hat{\theta}\) está 90% do tempo:

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

intervalo
amostras %>%
  mutate(no_intervalo = theta_c >= intervalo$valor_i &
           theta_c <= intervalo$valor_s) %>%
  summarise(cobertura = sum(no_intervalo) / n())
ggplot() +
  geom_rect(
    data = intervalo,
    aes(xmin = valor_i, xmax = valor_s),
    ymin = -Inf,
    ymax = Inf,
    fill = "coral", 
    alpha = .25
  ) +
  geom_histogram(
    data = amostras,
    aes(x = theta_c),
    binwidth = .5,
    fill = "white",
    colour = "darkgrey"
  ) +
  geom_vline(xintercept = THETA, linetype = "dashed") +
  labs(title = expression("Intervalo na distribuição amostral ao redor de" ~ theta))

Se sabemos a distância de A pra B, sabemos de B pra A

O pulo do gato é entender que quando sabemos dizer se \(\theta\) está perto de \(\hat{\theta}\), sabemos também dizer o inverso. Então se há 90% dos \(\hat{\theta}\) que estão entre \([\theta + \sigma_{.05}; \theta + \sigma_{.95}]\), isso também significa que para 90% dos \(\hat{\theta}\), \(\theta\) está a \([\hat{\theta} + \sigma_{.05}; \hat{\theta} + \sigma_{.95}]\).

Partindo dos \(\hat{\theta}\) que calculamos no início de nosso exercício, podemos ver como esse método funciona:

ic_amostras = amostras %>% 
  mutate(intervalo_i = theta_c + intervalo$erro_i, 
         intervalo_s = theta_c + intervalo$erro_s) %>% 
  mutate(contem_theta = THETA >= intervalo_i & THETA <= intervalo_s) 

ic_amostras %>% 
  sample_n(10)
ic_amostras %>% 
  summarise(cobertura = sum(contem_theta) / n())
ic_amostras %>% 
  sample_n(50) %>% 
  mutate(i = 1:n()) %>% 
  ggplot(aes(
    x = i,
    y = theta_c,
    ymin = intervalo_i,
    ymax = intervalo_s,
    color = contem_theta
  )) +
  geom_pointrange(alpha = .8, size = .3) +
  geom_hline(yintercept = THETA, color = "dark blue") +
  labs(x = "amostra",
       y = "média") +
  scale_color_manual(values = c("red", "grey70"))

Podia ser outra estatística

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

funcao_theta = function(df) {
  df %>%
    pull(news) %>%
    median()
}

THETA_OUTRO = funcao_theta(lastfm)

amostras = tibble(amostra = 1:1000) %>% # faremos 1000 vezes
  mutate(theta_c = map_dbl(amostra, ~ lastfm %>% 
                                       sample_n(100) %>%  
                                       funcao_theta()))

amostras
amostras %>% 
  ggplot(aes(theta_c)) + 
  geom_histogram(binwidth = 1, fill = "white", colour = "darkgrey") + 
  geom_vline(xintercept = THETA_OUTRO) 

amostras = amostras %>% 
  mutate(erro = theta_c - THETA_OUTRO)

amostras %>%
  ggplot(aes(erro)) +
  geom_histogram(binwidth = 1,
                 fill = "white",
                 colour = "darkblue") +
  geom_vline(xintercept = 0, linetype = "dashed") + 
  labs(title = "Distribuição do erro amostral")

margem_erro = amostras %>%
  summarise(baixo = quantile(erro, .05),
            cima = quantile(erro, .95)) 
margem_erro
amostras %>% 
  sample_n(1) %>% 
  mutate(ic_inferior = theta_c + margem_erro$baixo, 
         ic_superior = theta_c + margem_erro$cima)

Efeito do tamanho da amostra

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(artists) %>% 
                                      sample(100) %>% 
                                      mean()))

amostras %>% summary
##     amostra           media      
##  Min.   :   1.0   Min.   :18.72  
##  1st Qu.: 250.8   1st Qu.:21.30  
##  Median : 500.5   Median :22.08  
##  Mean   : 500.5   Mean   :22.13  
##  3rd Qu.: 750.2   3rd Qu.:22.91  
##  Max.   :1000.0   Max.   :26.29
amostras %>% 
  ggplot(aes(media)) + 
  geom_histogram(binwidth = 1, fill = "white", colour = "darkgrey") + 
  geom_vline(xintercept = THETA)

E como conseguimos a distribuição amostral?

Espero que até agora tenha ficado claro que tendo a distribuição dos valores de uma estatística \(\hat{\theta}\) a partir de amostras (a distribuição amostral) e o valor de \(\hat{\theta}\), nós conseguimos estimar um intervalo com um método que acerta com uma cobertura conhecida. Conseguimos estimar uma margem de erro a partir de nosso \(\hat{\theta}\) com uma confiança conhecida.

Só que…

Nós nunca temos várias amostras da população como nos exemplos até agora. Estamos aqui porque nós não temos a população, e todo o dado que temos forma uma amostra. Caso não tenha ficado claro: o exercício até agora foi simular o que aconteceria se tivéssemos várias amostras para fins pedagógicos.

A ideia central que usaremos

Não temos a população para estimar a distribuição amostral e a distribuição do erro amostral. Mas temos algo que veio dessa população, que é a amostra. Em várias situações como essa, a Estatística contorna a falta da informação ideal usando a que temos e um método que funcione bem assim.

A ideia principal que usaremos é uma técnica chamada boostrapping que funciona porque 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 que nos interessa (média, mediana, desvio padrão, o que for) para cada uma das \(b\) amostras, gerando \(b\) valores de \(\hat{\theta}*\).

Como resultado, teremos uma distribuição de como \(\hat{\theta}*\) varia.

O princípio do bootstrap diz que a variação de \(\hat{\theta}*\) nos bootstraps aproxima a variação de \(\hat{\theta}\).

Com isso, podemos usar a mesma lógica que usamos acima e construir um intervalo ao retor de \(\hat{\theta}\) que contém \(\theta\) com uma certa confiança.

Aplicando bootstrapping

\(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}]\).

Aplicando bootstrapping

glimpse(lastfm)
## Rows: 11,988
## Columns: 4
## $ news    <dbl> 38, 32, 24, 36, 28, 25, 30, 28, 14, 22, 27, 21, 23, 37, 25, 3…
## $ old     <dbl> 116, 130, 10, 46, 128, 95, 228, 53, 46, 118, 90, 34, 100, 109…
## $ ecletic <dbl> 1159.4960, 3439.8373, 1611.4464, 697.0974, 1660.1563, 1907.90…
## $ artists <dbl> 25.666667, 27.000000, 5.666667, 13.666667, 26.000000, 20.0000…
funcao_theta = function(df) {
  df %>%
    pull(artists) %>%
    mean()
}

THETA = funcao_theta(lastfm)

set.seed(1212)

amostra = lastfm %>%  
  sample_n(50) 

theta_c = funcao_theta(amostra)

theta_c
## [1] 24.08667
repeticoes = 4000 # pelo menos 2000, mas mais não faz mal.

um_bootstrap <- function(x){
  dados = x %>% pull(artists)
  boot_x <- sample(dados,              # amostre dos dados
                   size = NROW(dados), # tamanho igual ao recebido
                   replace = TRUE)     # aqui é o bootstrap
  return(mean(boot_x))
}

set.seed(1212)

# A REAMOSTRAGEM
reamostragens = tibble(i = 1:repeticoes) %>% 
  mutate(theta_c_s = map_dbl(i, ~ um_bootstrap(amostra)))

reamostragens
reamostragens %>%
  ggplot(aes(x = theta_c_s)) +
  geom_histogram(binwidth = .5,
                 colour = "darkorange",
                 fill = "white") + 
  labs(title = "Distribuição nos bootstraps", 
       subtitle = "É nossa estimativa da distribuição amostral")

reamostragens %>%
  ggplot(aes(x = theta_c_s - theta_c)) +
  geom_histogram(binwidth = 1,
                 colour = "darkblue",
                 fill = "white") + 
  labs(title = "Distribuição do erro amostral estimada via bootstraps")

Calculando o IC

Agora usamos a distribuição de \(\delta* = \hat{\theta}* - \hat{\theta}\) no lugar da de \(\delta\).

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

intervalo

Agora fazemos o mesmo que antes para estimar onde \(\theta\) está usando \(\hat{\theta}\).

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

intervalo
ggplot() +
  geom_rect(
    data = intervalo,
    aes(xmin = valor_i, xmax = valor_s),
    ymin = -Inf,
    ymax = Inf,
    fill = "gold",
    alpha = .25
  ) +
  annotate(
    geom = "text",
    x = intervalo$valor_i,
    y = -10,
    label = paste("hat(theta) + MoE_i"),
    parse = T,
    color = "brown",
    size = 3
  ) +
  annotate(
    geom = "text",
    x = intervalo$valor_s,
    y = -10,
    label = paste("hat(theta) + MoE_s"),
    parse = T,
    color = "brown",
    size = 3
  ) +
  geom_histogram(
    data = reamostragens,
    aes(theta_c_s),
    binwidth = .5,
    fill = "white",
    colour = "darkgrey"
  ) +
  geom_vline(xintercept = THETA,
             color = "dodgerblue",
             size = 1.2) +
  annotate(
    geom = "text",
    x = THETA - .3,
    y = -10,
    label = "theta",
    parse = T,
    color = "dodgerblue",
    size = 5
  ) +
  geom_vline(xintercept = theta_c, color = "brown") +
  annotate(
    geom = "text",
    x = theta_c - .3,
    y = -10,
    label = "hat(theta)",
    parse = T,
    color = "brown",
    size = 5
  ) +
  labs(title = expression("Intervalo estimado via bootstrap"))

Com outro nível de confiança:

confianca = .99
alpha = 1 - confianca

intervalo2 = reamostragens %>% 
  mutate(erro = theta_c_s - theta_c) %>% 
  summarise(erro_i = quantile(erro, alpha / 2), 
            erro_s = quantile(erro, 1 - alpha /2)) %>% 
  mutate(valor_i = theta_c + erro_i, 
         valor_s = theta_c + erro_s)

intervalo2
ggplot() +
  geom_rect(
    data = intervalo2,
    aes(xmin = valor_i, xmax = valor_s),
    ymin = -Inf,
    ymax = Inf,
    fill = "coral",
    alpha = .25
  ) +
  geom_rect(
    data = intervalo,
    aes(xmin = valor_i, xmax = valor_s),
    ymin = -Inf,
    ymax = Inf,
    fill = "coral2",
    alpha = .5
  ) +
  geom_histogram(
    data = reamostragens,
    aes(theta_c_s),
    binwidth = .5,
    fill = "white",
    colour = "darkgrey"
  ) +
  geom_vline(xintercept = THETA,
             color = "dodgerblue",
             size = 1.5) +
  geom_vline(xintercept = theta_c, color = "brown") +
  labs(title = expression("Intervalo estimado via bootstrap"), 
       subtitle = "Escuro: 99%, claro: 90% confiança", 
       x = expression(hat(theta) ~ "*"))

Como normalmente vemos esse IC

intervalo = intervalo %>% 
  mutate(theta_c = theta_c)

intervalo %>% 
  ggplot(aes(x = "", y = theta_c, ymin = valor_i, ymax = valor_s)) + 
  geom_linerange() + 
  geom_point(size = 2) + 
  coord_flip() +
  scale_y_continuous(limits = c(15, 30)) + 
  labs(
    title = "IC da média de artistas escutados", 
    y = "Média de artistas escutados por mês", 
    x = "")