Análise de uma replicação

Sobre IAT

IAT: 0.15, 0.35, and 0.65 are considered small, medium, and large level of bias for individual scores.

Positive means bias towards arts / against Math.

Análise de uma replicação

csv_files <- list.files(here::here("data"), pattern = "*.csv")
random_csv <- sample(csv_files, 1)

iats <- read_csv(here::here("data", random_csv), col_types = "cccdc") %>% 
  mutate(sex = factor(sex, levels = c("m", "f"), ordered = TRUE))

glimpse(iats)

## Rows: 125
## Columns: 5
## $ session_id  <chr> "2410996", "2411208", "2411213", "2413016", "2413017", "24…
## $ referrer    <chr> "ufl", "ufl", "ufl", "ufl", "ufl", "ufl", "ufl", "ufl", "u…
## $ sex         <ord> m, f, f, f, f, f, f, f, f, f, m, f, m, f, m, f, f, m, f, m…
## $ d_art       <dbl> 0.86629314, 0.68168619, 0.92370419, 1.23406708, 0.04755330…
## $ iat_exclude <chr> "Include", "Include", "Include", "Include", "Include", "In…

iats %>%
    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")

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

iats %>% 
    ggplot(aes(x = sex, y = d_art)) + 
    geom_quasirandom(width = .25) + 
    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 every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Qual a diferença na amostra

iats %>% 
  group_by(referrer, sex) %>% 
  summarise(media = round(mean(d_art), 3))

## `summarise()` has grouped output by 'referrer'. You can override using the
## `.groups` argument.

agrupado = iats %>% 
        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.3877554

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 = iats, 
               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.3877554
## $ bias      <dbl> 0.002206025
## $ std.error <dbl> 0.1106027
## $ conf.low  <dbl> -0.6093956
## $ conf.high <dbl> -0.1806341

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

p1 = iats %>% 
    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 = "Diferenca", 
         y = "IATS homens - mulheres")

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

Conclusão

Em média, as mulheres que participaram do experimento tiveram uma associação implícita (medida pelo IAT) com a matemárica positiva/negativa e fraca/média/forte (média 0.454, desv. padrão 0.12, N = 2000). Homens tiveram uma associação positiva/negativa com a matemática, portanto maior/menor que a das mulheres (média 0.081, desv. padrão 0.18, N = 2000). Houve portanto uma grande/considerável/pequena diferença entre homens e mulheres (diferença das médias -0.373, 95% CI \[-0.606, -0.149\]). A partir desta amostra, estimamos que 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 bibliotece (exemplo acima)

IC da média de IAT, na amostra de HOMENS Abington

library(boot)

theta_homens = function(d, i) {
    homens = iats %>%
        filter(sex == "m") %>%
        slice(i) %>%
        summarise(media_homens = mean(d_art))
    homens %>% pull(media_homens)
}

homens_booted <- boot(data = iats, 
                      statistic = theta_homens, 
                      R = 2000)

homens_ci = tidy(homens_booted, 
          conf.level = .95,
          conf.method = "bca",
          conf.int = TRUE)

homens_ci = cbind(homens_ci, sex = "m", desvio_padrao = (homens_ci$conf.high-homens_ci$conf.low)/2)
glimpse(homens_ci)

## Rows: 1
## Columns: 7
## $ statistic     <dbl> 0.03119136
## $ bias          <dbl> -0.000150545
## $ std.error     <dbl> 0.09471028
## $ conf.low      <dbl> -0.1492318
## $ conf.high     <dbl> 0.2318473
## $ sex           <chr> "m"
## $ desvio_padrao <dbl> 0.1905395

IC da média de IAT, na amostra de MULHERES Abington

library(boot)

theta_mulheres = function(d, i) {
    mulheres = iats %>%
        filter(sex == "f") %>%
        slice(i) %>%
        summarise(media_mulheres = mean(d_art))
    mulheres %>% pull(media_mulheres)
}

mulheres_booted <- boot(data = iats, 
                      statistic = theta_mulheres, 
                      R = 2000)

mulheres_ci = tidy(mulheres_booted, 
          conf.level = .95,
          conf.method = "bca",
          conf.int = TRUE)

mulheres_ci = cbind(mulheres_ci, sex = "f", desvio_padrao = (mulheres_ci$conf.high-mulheres_ci$conf.low)/2)
glimpse(mulheres_ci)

## Rows: 1
## Columns: 7
## $ statistic     <dbl> 0.4189468
## $ bias          <dbl> 0.0004106056
## $ std.error     <dbl> 0.05660811
## $ conf.low      <dbl> 0.3006261
## $ conf.high     <dbl> 0.5237714
## $ sex           <chr> "f"
## $ desvio_padrao <dbl> 0.1115726

Comparativo do IC da média de IAT, entre as amostras de HOMENS e MULHERES em Abington

mf_compar_ic = rbind(homens_ci, mulheres_ci)
mf_compar_ic %>%
ggplot(aes(x = sex, ymin = conf.low, y = statistic, ymax = conf.high)) +
  geom_errorbar(width = 0.2, position = position_dodge(0.5), color = "green") +
  geom_point(color = "green", size = 5) +
  geom_hline(yintercept = c(0.15, 0.35, 0.65), linetype = "dashed", color = "purple")

2. Bootstraps implementados por você (justifique o método de IC com bootstraps escolhido)

A função theta para calcular a média de uma amostra, agrupada pelo sexo:

funcao_theta <- function(df) {
  agrupado <- df %>%
    group_by(sex) %>% 
    summarise(media = mean(d_art))
  
  return(agrupado)
}

Obter a média da amostra, utilizando a função theta definida anteriormente, agrupada pelo sexo:

THETA <- funcao_theta(iats)
THETA_homens <- ( THETA %>% filter(sex == "m") %>% pull(media))
THETA_mulheres <- (THETA %>% filter(sex == "f") %>% pull(media))
THETA  # Media da amostra por sexo (homens x mulheres) 

A função de para gerar um bootstrap da media da amostra, agrupada pelo sexo:

my_bootstrap <- function(df) {
  dados_m <- df %>% filter(sex == "m") %>% pull(d_art)
  dados_f <- df %>% filter(sex == "f") %>% pull(d_art)
  boot_m <- sample(dados_m,
                   size = length(dados_m),
                   replace = TRUE)
  boot_f <- sample(dados_f,
                   size = length(dados_f),
                   replace = TRUE)
  return(c(media_m = mean(boot_m), media_f = mean(boot_f)))
}

Agora com 2000 que foi definido, gerando sub-amostras da media com bootstrap agrupada pelo sexo:

library(reshape2)

## 
## Attaching package: 'reshape2'

## The following object is masked from 'package:tidyr':
## 
##     smiths

amostras <- replicate(2000, my_bootstrap(iats), simplify = FALSE)
amostras <- do.call(rbind, amostras)
amostras <- as.data.frame(amostras)

amostras_melted <- melt(amostras, measure.vars = c("media_m", "media_f"), variable.name = "sexo")
amostras_melted$sexo <- factor(amostras_melted$sexo, labels = c("Homens", "Mulheres"))

densidade <- amostras_melted %>%
  ggplot(aes(x = value, fill = sexo)) +
  geom_density(alpha = 0.5) +
  labs(title = "Distribuicao bootstraps (Homens e Mulheres)",
       x = "Media") +
  scale_fill_manual(values = c("lightblue", "pink"))

densidade

hist_erro_m <- amostras %>%
  ggplot(aes(x = media_m - THETA_homens)) +
  geom_histogram(binwidth = .01,
                 colour = "black",
                 fill = "lightblue") + 
  labs(title = "Distribuição Erro Amostral (Homens)", x = "Diferença de Médias (M - THETA_homens)")

hist_erro_f <- amostras %>%
  ggplot(aes(x = media_f - THETA_mulheres)) +
  geom_histogram(binwidth = .01,
                 colour = "black",
                 fill = "pink") + 
  labs(title = "Distribuição Erro Amostral (Mulheres)", x = "Diferença de Médias (F - THETA_mulheres)")

grid.arrange(hist_erro_m, hist_erro_f, ncol = 2)

# Calcular a media para homens e mulheres separadamente
THETA_homens <- mean(iats %>% filter(sex == "m") %>% pull(d_art))
THETA_mulheres <- mean(iats %>% filter(sex == "f") %>% pull(d_art))

# Combinar os dados de homens e mulheres em um único data frame
dados <- rbind(
  data.frame(sex = "Homens", diff = amostras$media_m - THETA_homens),
  data.frame(sex = "Mulheres", diff = amostras$media_f - THETA_mulheres)
)

# Plotar o histograma do erro amostral para homens e mulheres no mesmo gráfico
hist_erro <- dados %>%
  ggplot(aes(x = diff, fill = sex)) +
  geom_histogram(binwidth = 0.01, colour = "black", alpha = 0.7) +
  scale_fill_manual(values = c("Homens" = "lightblue", "Mulheres" = "pink")) +
  labs(title = "Distribuicao Erro Amostral",
       x = "Diferenca de Medias",
       fill = "Sexo") +
  theme(legend.position = "top")

hist_erro

calcular o Intervalo de Confiança (IC) para homens e mulheres, com C = 95% e os limites inferior e superior da margem de erro (quantis: 2.5 (-) e 97.5 (+)).

Analisando a associação implícita entre matemática e gênero em Abington usando o IAT. No item anterior, aplicamos o método BCA com intervalo de confiança de 95%. Agora, reproduziremos os resultados com nosso método (my_bootstrap) para verificar mudanças na conclusão, mantendo o mesmo intervalo de confiança.

# Grafico de barras comparando as medias das associacoes implícitas entre homens e mulheres
barplot_data <- data.frame(sex = c("Homens", "Mulheres"),
                           statistic = c(THETA_homens, THETA_mulheres))

barplot <- ggplot(barplot_data, aes(x = sex, y = statistic, fill = sex)) +
  geom_bar(stat = "identity", width = 0.5) +
  labs(title = "Comparacao das medias das associacoes implicitas",
       x = "Genero", y = "Media") +
  scale_fill_manual(values = c("blue", "pink")) +
  theme_minimal()

# Exibir os histogramas e o grafico de barras juntos
grid.arrange(hist_erro_m, hist_erro_f, barplot, ncol = 2)

hist_m <- amostras %>%
  ggplot(aes(x = media_m)) +
  geom_histogram(binwidth = .01,
                 colour = "black",
                 fill = "lightblue") + 
  labs(title = "Distribuição Amostral (Homens)", x = "Média")

intervalo_m <- amostras %>% 
  mutate(erro_m = media_m - THETA_homens) %>% 
  summarise(erro_m_i = quantile(erro_m, .025), 
            erro_m_s = quantile(erro_m, .975)) %>%
  mutate(valor_m_i = THETA_homens + erro_m_i, 
         valor_m_s = THETA_homens + erro_m_s) %>%
  mutate(desvio_padrao_m = (valor_m_s - valor_m_i) / 2)

hist_comp_m <- hist_m +
  geom_vline(xintercept = c(intervalo_m$valor_m_i), linetype = "dashed", color = "darkblue") +
  geom_vline(xintercept = c(THETA_homens), size = 2, color = "red") +
  geom_vline(xintercept = c(intervalo_m$valor_m_s), linetype = "dashed", color = "darkgreen")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

hist_f <- amostras %>%
  ggplot(aes(x = media_f)) +
  geom_histogram(binwidth = .01,
                 colour = "black",
                 fill = "pink") + 
  labs(title = "Distribuição Amostral (Mulheres)", x = "Média")

intervalo_f <- amostras %>% 
  mutate(erro_f = media_f - THETA_mulheres) %>% 
  summarise(erro_f_i = quantile(erro_f, .025), 
            erro_f_s = quantile(erro_f, .975)) %>%
  mutate(valor_f_i = THETA_mulheres + erro_f_i, 
         valor_f_s = THETA_mulheres + erro_f_s) %>%
  mutate(desvio_padrao_f = (valor_f_s - valor_f_i) / 2)

hist_comp_f <- hist_f +
  geom_vline(xintercept = c(intervalo_f$valor_f_i), linetype = "dashed", color = "darkblue") +
  geom_vline(xintercept = c(THETA_mulheres), size = 2, color = "red") +
  geom_vline(xintercept = c(intervalo_f$valor_f_s), linetype = "dashed", color = "darkgreen")

grid.arrange(hist_comp_m, hist_comp_f, ncol = 2)

# Criar dataframe com intervalo de confiança para homens
intervalo_m_mycis <- data.frame(sex = "Homem", statistic = THETA_homens, conf.low = intervalo_m$valor_m_i, conf.high = intervalo_m$valor_m_s)

# Criar dataframe com intervalo de confiança para mulheres
intervalo_f_mycis <- data.frame(sex = "Mulher", statistic = THETA_mulheres, conf.low = intervalo_f$valor_f_i, conf.high = intervalo_f$valor_f_s)

# Combinar os dataframes de homens e mulheres
my_mf_cis <- rbind(intervalo_m_mycis, intervalo_f_mycis)

# Plotar o gráfico
my_mf_cis %>%
  ggplot(aes(x = sex, ymin = conf.low, y = statistic, ymax = conf.high)) +
  geom_linerange(linewidth = 2, color = "blue") +
  geom_point(color = "red", size = 5) +
  geom_hline(yintercept = c(0.15, 0.35, 0.65), linetype = "dashed", color = "darkgreen") +
  labs(x = "Sexo", y = "Estatistica", title = "Intervalo de Confianca") +
  scale_x_discrete(labels = c("Homem", "Mulher")) +  # Substituir os rotulos no eixo x
  theme_minimal()

Finalizar com a função bootstrap para calculando a diferença das médias

# Calcular a diferença das médias
diferenca_medias <- intervalo_m_mycis$statistic - intervalo_f_mycis$statistic

# Criar dataframe com a diferenca das médias e intervalo de confiança
diferenca_medias_mycis <- data.frame(sex = "Diferenca", statistic = diferenca_medias, conf.low = intervalo_m_mycis$conf.low - intervalo_f_mycis$conf.high, conf.high = intervalo_m_mycis$conf.high - intervalo_f_mycis$conf.low)

# Plotar o gráfico
diferenca_medias_mycis %>%
  ggplot(aes(x = sex, ymin = conf.low, y = statistic, ymax = conf.high)) +
  geom_linerange(linewidth = 2, color = "blue") +
  geom_point(color = "red", size = 5) +
  labs(x = "Sexo", y = "Diferenca das Medias", title = "Intervalo de Confianca para Diferenca das Medias") +
  scale_x_discrete(labels = c("Diferenca Media")) +  # Substituir o rotulo no eixo x
  theme_minimal()