ICs com o pacote boot

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 dois atributos que crio abaixo: o número (médio) de artistas novos (ie. previamente não escutados) escutados por mês pela pessoa, e se a pessoa é do grupo de 5% de pessoas que escuta artistas mais/menos populares.

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

amostra = lastfm %>% 
    filter(!is.na(news + old)) %>% 
    mutate(
        novos = news / 6, 
        grupo = case_when(
            mediana_pop >= quantile(mediana_pop, .95) ~ "fã de pop",
            mediana_pop <= quantile(mediana_pop, .05) ~ "não gosta de pop", 
            TRUE ~ "outros"
        )
    ) %>%
    filter(grupo != "outros") %>% 
    select(grupo, novos)

amostra %>% count(grupo)

glimpse(amostra)

## Rows: 1,200
## Columns: 2
## $ grupo <chr> "não gosta de pop", "não gosta de pop", "fã de pop", "fã de pop…
## $ novos <dbl> 4.166667, 4.666667, 3.500000, 3.833333, 2.000000, 20.166667, 23…

amostra %>% 
  ggplot(aes(x = novos)) + 
  facet_wrap(~ grupo, ncol = 1) + 
  geom_histogram(binwidth = 1, boundary = 0) + 
  geom_rug()

IC de uma estatística em um grupo

fas = amostra %>% 
  filter(grupo == "fã de pop")

fas %>% 
  summarise(media = mean(novos))

s <- function(d, i) {
    sumarizado = d %>% 
        slice(i) %>% 
        summarise(do_grupo = mean(novos))
    
    sumarizado %>% 
      pull(do_grupo)
}

s(fas, 1:(nrow(fas))) # theta_chapeu

## [1] 4.033333

library(boot)
library(broom)

booted <- boot(data = fas, 
               statistic = s, 
               R = 2000)

ci_fas = tidy(booted, 
              conf.level = .95,
              conf.method = "basic",
              conf.int = TRUE)

glimpse(ci_fas)

## Rows: 1
## Columns: 5
## $ statistic <dbl> 4.033333
## $ bias      <dbl> -0.0002558333
## $ std.error <dbl> 0.1071847
## $ conf.low  <dbl> 3.811688
## $ conf.high <dbl> 4.236375

ci_fas %>% 
  ggplot(aes(
            ymin = conf.low,
            y = statistic,
            ymax = conf.high,
            x = "Fãs de pop"
        )) +
        geom_linerange() +
        geom_point(color = "coral", size = 2) +
        scale_y_continuous(limits = c(0, 5)) +
        labs(x = "", y = "Artistas novos/mês") +
        coord_flip()

nao_fas = amostra %>% 
  filter(grupo != "fã de pop")

nao_fas %>% 
  summarise(media = mean(novos))

booted <- boot(data = nao_fas, 
               statistic = s, 
               R = 1000)
ci_naofas = tidy(booted, 
               conf.level = .95,
               conf.method = "basic",
               conf.int = TRUE)
glimpse(ci_naofas)

## Rows: 1
## Columns: 5
## $ statistic <dbl> 6.071389
## $ bias      <dbl> -0.009993889
## $ std.error <dbl> 0.2430288
## $ conf.low  <dbl> 5.573625
## $ conf.high <dbl> 6.518263

cis = bind_rows(
  "pop" = ci_fas, 
  "antipop" = ci_naofas, 
  .id = "grupo"
)

cis %>% 
  ggplot(aes(
            ymin = conf.low,
            y = statistic,
            ymax = conf.high,
            x = grupo
        )) +
        geom_linerange() +
        geom_point(color = "coral", size = 2) +
        scale_y_continuous(limits = c(2, 8)) +
        labs(x = "", y = "Artistas novos/mês") +
        coord_flip()

IC da diferença das médias

s <- function(d, i) {
    agrupado = d %>% 
        slice(i) %>% 
        group_by(grupo) %>% 
        summarise(do_grupo = mean(novos), .groups = "drop")
    a = agrupado %>% filter(grupo == "fã de pop") %>% pull(do_grupo)
    b = agrupado %>% filter(grupo == "não gosta de pop") %>% pull(do_grupo)
    a - b
}

theta_c = s(amostra, 1:nrow(amostra))

theta_c

## [1] -2.038056

booted <- boot(data = amostra, 
               statistic = s, 
               R = 2000)
ci = tidy(booted, 
          conf.level = .95,
          conf.method = "basic",
          conf.int = TRUE)
glimpse(ci)

## Rows: 1
## Columns: 5
## $ statistic <dbl> -2.038056
## $ bias      <dbl> -0.006059701
## $ std.error <dbl> 0.2594757
## $ conf.low  <dbl> -2.528261
## $ conf.high <dbl> -1.517547

ci %>% 
  ggplot(aes(
            ymin = conf.low,
            y = statistic,
            ymax = conf.high,
            x = ""
        )) +
        geom_linerange() +
        geom_point(color = "coral", size = 3) +
        scale_y_continuous(limits = c(-5, 5)) +
        labs(x = "", y = "Diferença na média (fãs - não fãs)") +
        coord_flip()