Exemplo de análise de uma replicação

iat = read_csv(here::here(params$arquivo_dados), col_types = "dccdc")
iat = iat %>% 
    mutate(sex = factor(sex, levels = c("m", "f"), ordered = TRUE))
glimpse(iat)

## Rows: 165
## Columns: 5
## $ session_id  <dbl> 2435230, 2435236, 2435237, 2435239, 2435240, 2435243, 2435…
## $ referrer    <chr> "swpson", "swpson", "swpson", "swpson", "swpson", "swpson"…
## $ sex         <ord> f, f, m, f, f, f, f, f, f, f, f, f, f, m, f, f, f, f, f, f…
## $ d_art       <dbl> 0.39171550, 1.00715226, 0.26248436, 0.64749897, 0.58174417…
## $ iat_exclude <chr> "Include", "Include", "Include", "Include", "Include", "In…

iat %>%
    ggplot(aes(x = d_art, fill = sex, color = sex)) +
    geom_histogram(binwidth = .2, alpha = .4) +
    geom_rug() +
    facet_grid(sex ~ ., scales = "free_y") + 
    theme(legend.position = "None")

iat %>% 
    ggplot(aes(x = sex, y = d_art)) + 
    geom_quasirandom(width = .1)

iat %>% 
    ggplot(aes(x = sex, y = d_art)) + 
    geom_quasirandom(width = .1) + 
    stat_summary(geom = "point", fun.y = "mean", color = "red", size = 5)

## Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
## ℹ Please use the `fun` argument instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Qual a diferença na amostra

iat %>% 
    group_by(sex) %>% 
    summarise(media = mean(d_art), dp = sd(d_art), n = n())

## # A tibble: 2 × 4
##   sex   media    dp     n
##   <ord> <dbl> <dbl> <int>
## 1 m     0.238 0.515    53
## 2 f     0.508 0.433   112

agrupado = iat %>% 
        group_by(sex) %>% 
        summarise(media = mean(d_art))
    m = agrupado %>% filter(sex == "m") %>% pull(media)
    f = agrupado %>% filter(sex == "f") %>% pull(media)
m - f

## [1] -0.2696448

Comparação via ICs

library(boot)

theta <- function(d, i) {
    agrupado = d %>% 
        slice(i) %>% 
        group_by(sex) %>% 
        summarise(media = mean(d_art))
    m = agrupado %>% filter(sex == "m") %>% pull(media)
    f = agrupado %>% filter(sex == "f") %>% pull(media)
    m - f
}

set.seed(1234)
booted <- boot(data = iat, 
               statistic = theta, 
               R = 5000)

ci = tidy(booted, 
          conf.level = .95,
          conf.method = "bca",
          conf.int = TRUE)

glimpse(ci)

## Rows: 1
## Columns: 5
## $ statistic <dbl> -0.2696448
## $ bias      <dbl> 0.0005327326
## $ std.error <dbl> 0.08168239
## $ conf.low  <dbl> -0.4322362
## $ conf.high <dbl> -0.1135285

ci %>%
    ggplot(aes(
        x = "",
        y = statistic,
        ymin = conf.low,
        ymax = conf.high
    )) +
    geom_pointrange() +
    geom_point(size = 3) + 
    labs(x = "Diferença", 
         y = "IAT homens - mulheres")

p1 = iat %>% 
    ggplot(aes(x = sex, y = d_art)) +
    geom_quasirandom(width = .1) + 
    stat_summary(geom = "point", fun.y = "mean", color = "red", size = 5)

p2 = ci %>%
    ggplot(aes(
        x = "",
        y = statistic,
        ymin = conf.low,
        ymax = conf.high
    )) +
    geom_pointrange() +
    geom_point(size = 3) + 
    ylim(-1, 1) + 
    labs(x = "Diferença", 
         y = "IAT homens - mulheres")

grid.arrange(p1, p2, ncol = 2)

Atividade

Preencha os resultados e conclusões abaixo

Em média, as mulheres que participaram do experimento tiveram uma associação implícita (medida pelo IAT) com a matemárica positiva e média (bem próxima do limiar com forte) (média 0.508, desv. padrão 0.433, N = 112). Homens tiveram uma associação positiva com a matemática, portanto menor que a das mulheres (média 0.238, desv. padrão 0.515, N = 53). Houve portanto uma diferença considerável entre homens e mulheres (diferença das médias -0.270, 95% CI [-0.432, -0.114]). A partir desta amostra, estimamos que a diferença real na população provavelmente está entre 0.11 e 0.43 ponto na escala IAT, faixa que cruza o limiar considerado pequeno (0.15) e chega perto do limiar médio (0.35), o suficiente para diferenciar uma diferença pequena de uma moderada, mas não para confirmar uma diferença grande.

Exemplos de possíveis conclusões para completar

mulheres têm uma associação negativa consideravelmente mais forte, com uma diferença que provavelmente está entre 0.6 e 1.0 ponto na escala IAT, o suficiente para diferenciar uma associação neutra de uma muito forte contra a matemática.
mulheres têm uma associação negativa mais forte, porém não é claro se essa diferença é grande, moderada ou pequena. É necessário coletar mais dados para determinar se a diferença é relevante ou negligenciável.
mulheres podem ter um associação negativa forte, pode não haver diferença entre sexos ou homens podem ter atitudes negativas um pouco mais fortes ou moderadamente mais fortes.
pode não haver uma diferença entre sexos, ou se ela existir, ela provavelmente é pequena em qualquer das direções.

Realize novas análises sobre IAT usando as abordagens a seguir

Realize a análise e compare as conclusões obtidas nos dois casos experimentados:

bootstraps a partir de uma biblioteca (exemplo acima)
bootstraps implementados por você utilizando pelo menos dois outros métodos de IC com bootstrap
justifique e escolha dos métodos e discuta / explique caso exista impacto nos resultados

1. Checando exclusões do estudo original

Antes de reamostrar, verificamos se o critério de exclusão do estudo original marcou algum caso como inválido.

count(iat, iat_exclude)

## # A tibble: 1 × 2
##   iat_exclude     n
##   <chr>       <int>
## 1 Include       165

Todos os 165 casos estão marcados como "Include". Filtramos explicitamente por esse valor para que o código permaneça correto caso a base mude:

iat = iat %>%
    filter(iat_exclude == "Include")

nrow(iat)

## [1] 165

2. Bootstrap implementado manualmente

Implementamos o próprio algoritmo de reamostragem, sem usar boot(): cada reamostra seleciona 165 observações com reposição (mesmo tamanho da amostra original), recalcula a diferença de médias e armazena o resultado.

bootstrap_manual <- function(dados, n_boot = 5000, seed = 1234) {
    set.seed(seed)
    n <- nrow(dados)

    map_dfr(1:n_boot, function(i) {
        amostra <- dados %>% slice_sample(n = n, replace = TRUE)

        medias <- amostra %>%
            group_by(sex) %>%
            summarise(media = mean(d_art), .groups = "drop")

        m <- medias %>% filter(sex == "m") %>% pull(media)
        f <- medias %>% filter(sex == "f") %>% pull(media)

        tibble(media_m = m, media_f = f, diferenca = m - f)
    })
}

boot_manual <- bootstrap_manual(iat)

boot_manual %>%
    ggplot(aes(x = diferenca)) +
    geom_histogram(bins = 40, fill = "steelblue", alpha = .7) +
    geom_vline(xintercept = mean(boot_manual$diferenca), color = "red") +
    geom_vline(xintercept = m - f, color = "darkgreen", linetype = "dashed") +
    labs(
        title = "Distribuição bootstrap da diferença (m - f)",
        subtitle = "linha vermelha: média do bootstrap | linha verde: estimativa na amostra original",
        x = "diferença de médias (bootstrap)",
        y = "contagem"
    )

A distribuição acima é aproximadamente simétrica, com uma média do bootstrap (-0.2692) muito próxima da estimativa na amostra original (-0.2696), ligeiramente menos negativa. Essa pequena discrepância é justamente o viés que os diferentes métodos de IC tratam de formas diferentes, o que motiva comparar mais de um método em vez de usar só um.

A partir dessa mesma distribuição bootstrap, calculamos os intervalos por dois métodos distintos:

Método A: Percentil

O método mais direto: os limites do IC são os percentis 2.5% e 97.5% da própria distribuição bootstrap. Assume que a distribuição bootstrap já aproxima bem a distribuição amostral da estatística, sem correções adicionais.

ic_percentil <- boot_manual %>%
    summarise(
        estimate  = mean(diferenca),
        conf.low  = quantile(diferenca, .025),
        conf.high = quantile(diferenca, .975)
    )

ic_percentil

## # A tibble: 1 × 3
##   estimate conf.low conf.high
##      <dbl>    <dbl>     <dbl>
## 1   -0.269   -0.427    -0.112

Método B: Basic (reflexão em torno da estimativa observada)

Diferente do percentil, o método Basic não usa os percentis da distribuição bootstrap diretamente como limites do IC. Em vez disso, reflete esses percentis em torno da estimativa observada na amostra original (-0.2696), com a lógica de que o desvio entre a reamostra e a estimativa original (d* = θ* − θ̂) é o que deveria ser usado para construir o intervalo, não o valor absoluto de θ*:

\[IC_{basic} = \left[\,2\hat\theta - q_{0.975},\;\; 2\hat\theta - q_{0.025}\,\right]\]

onde q_{0.025} e q_{0.975} são os mesmos percentis da distribuição bootstrap usados no método Percentil.

theta_hat <- m - f

ic_basic <- tibble(
    estimate  = theta_hat,
    conf.low  = 2 * theta_hat - quantile(boot_manual$diferenca, .975),
    conf.high = 2 * theta_hat - quantile(boot_manual$diferenca, .025)
)

ic_basic

## # A tibble: 1 × 3
##   estimate conf.low conf.high
##      <dbl>    <dbl>     <dbl>
## 1   -0.270   -0.428    -0.112

3. Comparação entre os três métodos

comparacao <- bind_rows(
    ci %>%
        transmute(metodo = "Biblioteca (BCa)",
                  estimate = statistic, conf.low, conf.high),
    ic_percentil %>% mutate(metodo = "Manual (Percentil)"),
    ic_basic %>% mutate(metodo = "Manual (Basic)")
) %>%
    select(metodo, estimate, conf.low, conf.high) %>%
    mutate(largura_ic = conf.high - conf.low)

comparacao

## # A tibble: 3 × 5
##   metodo             estimate conf.low conf.high largura_ic
##   <chr>                 <dbl>    <dbl>     <dbl>      <dbl>
## 1 Biblioteca (BCa)     -0.270   -0.432    -0.114      0.319
## 2 Manual (Percentil)   -0.269   -0.427    -0.112      0.316
## 3 Manual (Basic)       -0.270   -0.428    -0.112      0.316

comparacao %>%
    ggplot(aes(x = metodo, y = estimate, ymin = conf.low, ymax = conf.high)) +
    geom_pointrange(size = 1) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "grey40") +
    coord_flip() +
    labs(
        title = "IC 95% para a diferença (homens - mulheres): três métodos",
        x = NULL,
        y = "diferença de médias no IAT"
    )

Discussão da comparação

Com os dados reais (BCa: -0.270, IC [-0.432, -0.114]; Percentil: -0.269, IC [-0.427, -0.112]; Basic: -0.270, IC [-0.428, -0.112]), os três métodos concordam na direção (todos negativos, nenhum cruza o zero) e têm larguras de intervalo muito próximas (0.319 para o BCa contra 0.316 para os dois métodos manuais), porque o tamanho amostral (N = 165) já é grande o suficiente para que o viés e a assimetria da distribuição bootstrap sejam modestos.

A diferença que existe entre eles não é de largura, e sim de posicionamento: o método Basic desloca o intervalo em relação ao Percentil na direção oposta ao viés observado na distribuição bootstrap. Como a média do bootstrap ficou ligeiramente menos negativa do que a estimativa da amostra original, o método Basic corrige isso deslocando o intervalo para valores um pouco mais negativos, enquanto o Percentil aceita os percentis da distribuição bootstrap sem essa correção. O BCa da biblioteca faz uma correção semelhante à do Basic, mas soma a ela uma correção adicional de assimetria (o parâmetro de aceleração, estimado por jackknife), o que explica por que o intervalo do BCa é um pouco mais largo que o do Basic, mesmo corrigindo o mesmo tipo de viés.

Dado que a distribuição bootstrap desta amostra mostra apenas uma assimetria leve, os três métodos convergem para conclusões práticas equivalentes aqui. Isso não seria necessariamente verdade em uma amostra menor ou com uma distribuição mais assimétrica, onde a escolha do método poderia mudar a conclusão sobre se o IC cruza ou não um limiar relevante (por exemplo, o zero, ou os limiares de efeito pequeno/médio/grande do IAT). Diante disso, o BCa é o método que inspira mais confiança para generalizar esta estimativa, por corrigir tanto o viés quanto a assimetria observados; os métodos Percentil e Basic servem aqui como verificação de robustez, e o fato de os três concordarem é evidência de que a conclusão não depende de detalhes de implementação.

Análise de uma replicação

Sobre Implicit Association Tests (IAT) – Olhe o README do repositório.