Análise de uma replicação

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

IAT: 0.15, 0.35, and 0.65 are considered small, medium, and large levels of bias for individual scores. Positive means bias towards arts / against Math.

Os dados usados aqui vêm de uma replicação, conduzida em cerca de 30 laboratórios ao redor do mundo, do estudo de Nosek, Banaji & Greenwald (2002), “Sex differences in implicit math attitudes”. Os dados tratados foram disponibilizados por Robert J. Calin-Jageman para o workshop Teaching the New Statistics (osf.io/wx7ck), a partir dos dados originais da replicação (osf.io/pqf9r).

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:

1. bootstraps a partir de uma biblioteca (exemplo acima)
2. bootstraps implementados por você utilizando pelo menos dois outros métodos de IC com bootstrap
3. 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 r nrow(iat) casos estão marcados como “Include”. Foi aplicado um filtro explícito 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

m, f e a estimativa pontual theta_hat foram recalculados depois do filtro acima, e não a partir dos valores calculados na seção anterior. Nesta base específica não há exclusões, então isso não muda os números; entretanto, se a base usada tivesse casos marcados como “Exclude”, o uso do m/f antigos (calculados antes do filtro) faria com que a estimativa pontual usada no IC não correspondesse aos dados de fato reamostrados. O cálculo a partir do iat já filtrado deixa o código correto independentemente da base:

agrupado = iat %>%
    group_by(sex) %>%
    summarise(media = mean(d_art), .groups = "drop")

m = agrupado %>% filter(sex == "m") %>% pull(media)
f = agrupado %>% filter(sex == "f") %>% pull(media)
theta_hat = m - f

theta_hat

## [1] -0.2696448

2. Bootstrap implementado manualmente

Foi impelementado o próprio algoritmo de reamostragem, sem usar boot(). Feito de duas formas, para checar se o desenho da reamostragem importa aqui:

• Não-estratificado: cada reamostra seleciona 165 observações com reposição do conjunto inteiro (mesmo tamanho da amostra original), do mesmo jeito que o exemplo com boot() acima faz implicitamente. Isso significa que o número de homens e mulheres pode variar de reamostra para reamostra (nesta amostra original há N=53 homens e N=112 mulheres, mas uma reamostra poderia sair, por exemplo, com 48 homens e 117 mulheres).
• Estratificado por grupo: cada reamostra reamostra homens e mulheres separadamente, sempre preservando os tamanhos originais dos grupos (53 homens, 112 mulheres). Essa é a forma mais comum de se fazer bootstrap ao comparar dois grupos independentes, já que trata o tamanho de cada grupo como fixo (parte do desenho do estudo), e não como algo que deveria variar aleatoriamente.

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

    map_dfr(1:n_boot, function(i) {
        if (estratificado) {
            amostra <- dados %>%
                group_by(sex) %>%
                slice_sample(prop = 1, replace = TRUE) %>%
                ungroup()
        } else {
            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, estratificado = FALSE)
boot_estratificado <- bootstrap_manual(iat, estratificado = TRUE)

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 = theta_hat, 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.

Checagem: nossa implementação bate com a do boot()?

Antes do cálculo dos ICs, a distribuição gerada pela função manual foi comparada com a distribuição gerada pelo boot() da biblioteca (booted$t), para confirmar que as duas implementações geram, de fato, distribuições equivalentes, e não apenas produzem ICs parecidos por coincidência:

bind_rows(
    tibble(diferenca = as.numeric(booted$t), origem = "boot() (biblioteca)"),
    tibble(diferenca = boot_manual$diferenca, origem = "bootstrap_manual() (não-estratificado)")
) %>%
    ggplot(aes(x = diferenca, fill = origem)) +
    geom_density(alpha = .4) +
    labs(
        title = "Distribuição bootstrap: biblioteca vs. implementação manual",
        x = "diferença de médias (bootstrap)",
        y = "densidade",
        fill = NULL
    ) +
    theme(legend.position = "bottom")

tibble(
    origem = c("boot() (biblioteca)", "bootstrap_manual() (não-estratificado)"),
    media  = c(mean(booted$t), mean(boot_manual$diferenca)),
    dp     = c(sd(booted$t), sd(boot_manual$diferenca))
)

## # A tibble: 2 × 3
##   origem                                  media     dp
##   <chr>                                   <dbl>  <dbl>
## 1 boot() (biblioteca)                    -0.269 0.0817
## 2 bootstrap_manual() (não-estratificado) -0.269 0.0816

As duas distribuições se sobrepõem quase completamente, com média e desvio padrão muito próximos entre si. Isso é esperado, já que as duas usam o mesmo desenho de reamostragem (reposição sobre o conjunto inteiro, sem estratificação por grupo), sendo a semente e o gerador de números aleatórios usados internamente por cada implementação a única diferença entre elas. Essa checagem confirma que bootstrap_manual() está correta, antes de seu uso para o cálculo dos ICs abaixo.

A partir da distribuição bootstrap não-estratificada, os intervalos foram calculados 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.

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

Método C: Bootstrap-t (studentizado)

Percentil e Basic, nesta amostra, tendem a concordar bastante porque a distribuição bootstrap é aproximadamente simétrica, sendo a diferença entre eles apenas de posicionamento, não de largura. Para disponibilizar um método com maior chance de discordar de fato dos outros dois (permitindo testar se a conclusão realmente independe do método escolhido, e não apenas desta amostra em particular), o bootstrap-t também foi implementado.

A ideia do bootstrap-t é “estudentizar” cada reamostra: em vez de olhar para os percentis de θ* diretamente, olha-se para os percentis da estatística t* = (θ* − θ̂) / SE, onde SE é o erro padrão da própria reamostra. Isso exige um SE por reamostra; por isso, foi utilizado um bootstrap aninhado (bootstrap dentro do bootstrap) para estimar o SE* de cada reamostra externa:

se_diferenca <- function(dados) {
    dados %>%
        group_by(sex) %>%
        summarise(media = mean(d_art), variancia = var(d_art), n = n(),
                  .groups = "drop") %>%
        summarise(se = sqrt(sum(variancia / n))) %>%
        pull(se)
}

se_hat <- se_diferenca(iat)

bootstrap_t <- function(dados, n_boot = 2000, n_boot_interno = 200, 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_estrela <- medias %>% filter(sex == "m") %>% pull(media)
        f_estrela <- medias %>% filter(sex == "f") %>% pull(media)
        diferenca_estrela <- m_estrela - f_estrela

        se_estrela <- se_diferenca(amostra)

        tibble(diferenca = diferenca_estrela,
               t_estrela = (diferenca_estrela - theta_hat) / se_estrela)
    })
}

boot_t <- bootstrap_t(iat)

q_t <- quantile(boot_t$t_estrela, c(.025, .975), na.rm = TRUE)

ic_boot_t <- tibble(
    estimate  = theta_hat,
    conf.low  = theta_hat - q_t[2] * se_hat,
    conf.high = theta_hat - q_t[1] * se_hat
)

ic_boot_t

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

Diferente do Percentil e do Basic, o bootstrap-t não utiliza os quantis de θ* na escala original: os quantis utilizados são os da estatística estudentizada t, usada para reconstruir o intervalo na escala de θ̂. Em teoria, isso o torna mais preciso (assintoticamente mais próximo do nível de confiança nominal) que Percentil e Basic, especialmente quando a distribuição amostral de θ̂ não é bem aproximada por uma normal, ao custo de ser computacionalmente mais caro (por causa do bootstrap aninhado) e de poder gerar ICs mais instáveis quando SE fica muito próximo de zero em alguma reamostra. Nesta amostra, o IC obtido foi de -0.270, IC 95% [-0.436, -0.116], intervalo ligeiramente mais largo (0.320) que o dos demais métodos manuais (0.316) e até um pouco mais largo que o do BCa (0.319).

3. Comparação entre os quatro 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)"),
    ic_boot_t %>% mutate(metodo = "Manual (Bootstrap-t)")
) %>%
    select(metodo, estimate, conf.low, conf.high) %>%
    mutate(largura_ic = conf.high - conf.low)

comparacao

## # A tibble: 4 × 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
## 4 Manual (Bootstrap-t)   -0.270   -0.436    -0.116      0.320

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): quatro métodos",
        x = NULL,
        y = "diferença de médias no IAT"
    )

Adicionalmente, foi comparado a estimativa pontual e a largura do IC (via Percentil) obtidas com a reamostragem não-estratificada contra a reamostragem estratificada por grupo, para checar se o desenho da reamostragem (seção 2) muda a conclusão:

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

bind_rows(
    ic_percentil %>% mutate(metodo = "Percentil (não-estratificado)"),
    ic_percentil_estratificado %>% mutate(metodo = "Percentil (estratificado por grupo)")
) %>%
    mutate(largura_ic = conf.high - conf.low) %>%
    select(metodo, estimate, conf.low, conf.high, largura_ic)

## # A tibble: 2 × 5
##   metodo                              estimate conf.low conf.high largura_ic
##   <chr>                                  <dbl>    <dbl>     <dbl>      <dbl>
## 1 Percentil (não-estratificado)         -0.269   -0.427    -0.112      0.316
## 2 Percentil (estratificado por grupo)   -0.268   -0.432    -0.106      0.326

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]; Bootstrap-t: -0.270, IC [-0.436, -0.116]), os quatro métodos concordam na direção (todos negativos, nenhum cruza o zero) e apresentam larguras de intervalo próximas entre si (0.319 para o BCa, 0.316 para Percentil e Basic, e 0.320 para o Bootstrap-t), pois 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 entre Percentil e Basic não é de largura, mas de posicionamento: o intervalo do método Basic é deslocado, 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 que a estimativa da amostra original, essa diferença é corrigida pelo método Basic por meio de um deslocamento do intervalo para valores um pouco mais negativos, enquanto no Percentil os percentis da distribuição bootstrap são aceitos sem essa correção. Uma correção semelhante à do Basic é feita pelo BCa da biblioteca, à qual se soma 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.

O método com maior potencial de discordância em relação aos demais é o bootstrap-t, já que o intervalo é reconstruído a partir da distribuição da estatística estudentizada t, e não dos quantis de θ diretamente. Nesta amostra, o intervalo do bootstrap-t (largura 0.320) foi de fato o mais largo entre os quatro métodos, ainda que por uma margem pequena, o que é coerente com a maior sensibilidade desse método a reamostras com erro padrão pequeno. A conclusão prática, no entanto, não é alterada: mesmo o limite superior mais próximo de zero, obtido pelo bootstrap-t (-0.116), permanece bem distante de cruzar o zero.

A comparação entre os desenhos de reamostragem também não altera a conclusão. O IC obtido pelo Percentil com reamostragem estratificada por grupo (-0.268, [-0.432, -0.106], largura 0.326) fica muito próximo do obtido com a reamostragem não-estratificada (-0.269, [-0.427, -0.112], largura 0.316), com uma diferença de largura pequena (cerca de 3%) e limites superiores e inferiores deslocados por menos de 0.005, respectivamente. A leve ampliação do intervalo no desenho estratificado é coerente com o esperado: ao fixar os tamanhos dos grupos (53 homens, 112 mulheres) em vez de permitir que variem entre reamostras, uma fonte adicional de variação é removida da reamostragem simples, mas essa fonte de variação tende a ser pequena quando os grupos já têm tamanho razoável, como é o caso aqui. Em ambos os desenhos, o zero permanece fora do intervalo por uma margem confortável, reforçando que a conclusão sobre a diferença entre os sexos não depende de qual dos dois esquemas de reamostragem é adotado.

De modo geral, os quatro métodos convergem para conclusões práticas equivalentes nesta amostra, o que é esperado dado que a distribuição bootstrap apresenta apenas assimetria leve. Em uma amostra menor ou com distribuição mais assimétrica, o mesmo não poderia ser assumido: a escolha do método e do desenho da reamostragem poderia alterar 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). O método que inspira mais confiança para a generalização desta estimativa continua sendo o BCa, por corrigir tanto o viés quanto a assimetria observados; os métodos Percentil, Basic e Bootstrap-t funcionam aqui como verificação de robustez, e a concordância entre os quatro, somada à checagem realizada contra a implementação da biblioteca, é evidência de que a conclusão não depende de detalhes de implementação nem do método de IC escolhido.