Análise de uma replicação

iat = read_csv(here::here(params$arquivo_dados), col_types = "cccdc")
iat = iat %>% 
    mutate(sex = factor(sex, levels = c("m", "f"), ordered = TRUE))
glimpse(iat)
## Rows: 155
## Columns: 5
## $ session_id  <chr> "2436706", "2436967", "2440429", "2440430", "2440431", "24…
## $ referrer    <chr> "sdsu", "sdsu", "sdsu", "sdsu", "sdsu", "sdsu", "sdsu", "s…
## $ sex         <ord> f, f, f, f, m, f, f, m, f, m, f, f, f, f, f, f, m, m, f, m…
## $ d_art       <dbl> 0.90444320, -0.47402625, 0.46840862, -0.02522412, 0.136813…
## $ 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) + 
    stat_summary(geom = "point", fun.y = "mean", color = "red")
## 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 every 8 hours.
## 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),
              desv_padrao = sd(d_art),
              N = length(d_art))
## # A tibble: 2 × 4
##   sex   media desv_padrao     N
##   <ord> <dbl>       <dbl> <int>
## 1 m     0.224       0.485    38
## 2 f     0.467       0.548   117
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.2430539

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
}

booted <- boot(data = iat, 
               statistic = theta, 
               R = 2000)

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

glimpse(ci)
## Rows: 1
## Columns: 5
## $ statistic <dbl> -0.2430539
## $ bias      <dbl> -0.001337485
## $ std.error <dbl> 0.09342751
## $ conf.low  <dbl> -0.4319103
## $ conf.high <dbl> -0.06613866
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)

Conclusão

Em média, as mulheres que participaram do experimento apresentaram uma tendência, medida pelo IAT, positiva e moderada às artes (média = 0.4666898, desvio padrão = 0.5475448 e N = 117), cabe destacar o alto desvio padrão que indica que talvez a média não represente bem esta amostra. Diferentemente nos homens, em média, não parece haver preferência por nenhuma das áreas (média = 0.2236359, desvio padrão = 0.5475448, N = 38, observamos novamente um alto desvio padrão. Houve portanto uma diferença entre mulheres e homens (diferença de médias de -0.243, intervalo de confiança de 95% [-0.425, -0.056]). No entanto, a magnitude dessa diferença não pode ser determinada com exatidão, pois o intervalo de confiança abrange valores considerados moderados e elevados. Assim, os dados de nosso experimento apontam que mulheres tem uma maior associação com as artes que os homens, 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.

1. bootstraps a partir de uma biblioteca

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
}

booted <- boot(data = iat, 
               statistic = theta, 
               R = 10000)

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

print(ci)
## # A tibble: 1 × 5
##   statistic     bias std.error conf.low conf.high
##       <dbl>    <dbl>     <dbl>    <dbl>     <dbl>
## 1    -0.243 -0.00215    0.0928   -0.425   -0.0643

2. bootstraps manual

# Função para calcular a estatística de interesse (diferença de médias)
theta <- function(data, indices) {
  amostra <- data[indices, ]
  agrupado <- aggregate(d_art ~ sex, data = amostra, FUN = mean)
  m <- agrupado[agrupado$sex == "m", "d_art"]
  f <- agrupado[agrupado$sex == "f", "d_art"]
  m - f
}

# Número de amostras Bootstrap
n_amostras <- 10000

# Vetor para armazenar as estatísticas Bootstrap
estatisticas_bootstrap <- numeric(n_amostras)

# Loop para gerar as amostras Bootstrap e calcular a estatística de interesse
set.seed(123)  # Define a semente para garantir a reprodutibilidade
for (i in 1:n_amostras) {
  amostra_indices <- sample(nrow(iat), replace = TRUE)
  estatisticas_bootstrap[i] <- theta(iat, amostra_indices)
}

# Ordenar as estatísticas Bootstrap
estatisticas_bootstrap_ord <- sort(estatisticas_bootstrap)

# Intervalo de confiança Bootstrap de 95% (percentis)
percentil_inferior <- 0.025
percentil_superior <- 0.975

posicao_inferior <- floor(n_amostras * percentil_inferior) + 1
posicao_superior <- floor(n_amostras * percentil_superior) + 1

limite_inferior <- estatisticas_bootstrap_ord[posicao_inferior]
limite_superior <- estatisticas_bootstrap_ord[posicao_superior]

# Imprimir intervalo de confiança Bootstrap
intervalo_confianca <- c(limite_inferior, limite_superior)

# Cálculo do desvio padrão
desvio_padrao <- sd(estatisticas_bootstrap)

print(intervalo_confianca)
## [1] -0.42739909 -0.06181045
print(desvio_padrao)
## [1] 0.09331206

Resultados

Nesses ultimos casos o método de Intervalo de Confiança (IC) escolhido foi o método “perc”, que se refere ao Bootstrap Percentile. Esse método é especialmente útil quando os dados têm uma distribuição assimétrica ou quando os estimadores são enviesados. Ele utiliza percentis da distribuição empírica para capturar essa assimetria, fornecendo intervalos de confiança robustos. Além disso, o número de de amostras foi incrementado para 10000, afim de analisar mais casos. Contudo, temos praticamente o mesmo resultado do mostrado no exemplo com 4000 amostras, onde as que mulheres tem uma maior associação com as artes que os homens, 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.