Lab 04 - FPCC 2

Para o presente laboratório foi escolhido a seguinte fonte de dados: data/vcu .csv. A partir dos dados, é possível observar as sguintes visualizações:

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(ggplot2)
library(broom)
library(ggbeeswarm)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
library(boot)

theme_set(theme_bw())

knitr::opts_chunk$set(tidy = FALSE,
                      fig.width = 6,
                      fig.height = 5)

dados = read_csv(here::here("data/vcu .csv"),)
## Rows: 101 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): referrer, sex, iat_exclude
## dbl (2): session_id, d_art
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
glimpse(dados)
## Rows: 101
## Columns: 5
## $ session_id  <dbl> 2448165, 2451903, 2451987, 2457829, 2458789, 2458795, 2458…
## $ referrer    <chr> "vcu", "vcu", "vcu", "vcu", "vcu", "vcu", "vcu", "vcu", "v…
## $ sex         <chr> "f", "f", "f", "m", "f", "f", "f", "m", "m", "f", "f", "f"…
## $ d_art       <dbl> 1.074677138, 0.293409771, 0.006185518, -0.116358515, 1.011…
## $ iat_exclude <chr> "Include", "Include", "Include", "Include", "Include", "In…
iat = dados %>% 
    mutate(sex = factor(sex, levels = c("m", "f"), ordered = TRUE))
glimpse(iat)
## Rows: 101
## Columns: 5
## $ session_id  <dbl> 2448165, 2451903, 2451987, 2457829, 2458789, 2458795, 2458…
## $ referrer    <chr> "vcu", "vcu", "vcu", "vcu", "vcu", "vcu", "vcu", "vcu", "v…
## $ sex         <ord> f, f, f, m, f, f, f, m, m, f, f, f, f, f, f, f, m, f, f, f…
## $ d_art       <dbl> 1.074677138, 0.293409771, 0.006185518, -0.116358515, 1.011…
## $ 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 every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

iat %>% 
    group_by(sex) %>% 
    summarise(media = mean(d_art))
## # A tibble: 2 × 2
##   sex    media
##   <ord>  <dbl>
## 1 m     0.0664
## 2 f     0.468
agrupado = iat %>% 
    group_by(sex) %>% 
    summarise(media = mean(d_art))
m = agrupado %>% filter(sex == "m") %>% pull(media)
f = agrupado %>% filter(sex == "f") %>% pull(media)
f - m
## [1] 0.401629
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)
    f - m
}

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.401629
## $ bias      <dbl> 0.0006519254
## $ std.error <dbl> 0.1259744
## $ conf.low  <dbl> 0.1475587
## $ conf.high <dbl> 0.6481789
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")

plot_with_statistics <- function(data) {
    # Calculando as estatísticas
    mean_value <- mean(data$d_art, na.rm = TRUE)
    sd_value <- sd(data$d_art, na.rm = TRUE)
    
    # Calculando o intervalo de confiança de 95%
    n <- sum(!is.na(data$d_art))
    error_margin <- qnorm(0.975) * (sd_value / sqrt(n))
    ci_lower <- mean_value - error_margin
    ci_upper <- mean_value + error_margin
    
    # Criando um data frame para o ggplot
    stats_df <- data.frame(
        Statistic = "Statistics",
        Mean = mean_value,
        SD = sd_value,
        CI_Lower = ci_lower,
        CI_Upper = ci_upper
    )

    # Gerando o gráfico
    dpr <-  ggplot(stats_df, aes(x = Statistic, y = Mean)) +
        geom_col(fill = "#69b3a2", width = 0.5) +  # Usando um tom pastel
        geom_errorbar(aes(ymin = CI_Lower, ymax = CI_Upper), width = 0.2, color = "#404080") +  # Tons pastéis
        labs(title = "Mean with 95% Confidence Interval",
             x = "",
             y = "Value") +
        theme_minimal() +
        theme(
            plot.title = element_text(size = 15, face = "bold", hjust = 0.5),
            axis.title.y = element_text(size = 12),
            axis.text = element_text(size = 10)
        ) +
        geom_text(aes(label = round(Mean, 2)), vjust = -1.5, hjust = 1.5, color = "#404080", size = 5) +  # Adicionando valor da média
        geom_point(aes(y = Mean), size = 4, color = "#404080")  # Adicionando ponto da média
    
    return(dpr)
}

plot_with_statistics_by_gender <- function(data, column_name, gender_column) {
    # Filtrando e agrupando os dados por gênero
    stats_df <- data %>%
        group_by(!!sym(gender_column)) %>%
        summarise(
            Mean = mean(!!sym(column_name), na.rm = TRUE),
            SD = sd(!!sym(column_name), na.rm = TRUE),
            N = sum(!is.na(!!sym(column_name)))
        ) %>%
        mutate(
            Error_Margin = qnorm(0.975) * (SD / sqrt(N)),
            CI_Lower = Mean - Error_Margin,
            CI_Upper = Mean + Error_Margin
        )
    
    # Determinando os limites do eixo y, incluindo possíveis valores negativos
    y_min <- min(stats_df$CI_Lower, na.rm = TRUE)
    y_max <- max(stats_df$CI_Upper, na.rm = TRUE)
    
    # Gerando o gráfico
    dpr <- ggplot(stats_df, aes(x = !!sym(gender_column), y = Mean, fill = !!sym(gender_column))) +
        geom_col(width = 0.5) +
        geom_errorbar(aes(ymin = CI_Lower, ymax = CI_Upper), width = 0.2, color = "gray30") +
        scale_fill_manual(values = c("m" = "#69b3a2", "f" = "#ff9999")) +  # Usando tons pastéis
        labs(title = "Mean with 95% Confidence Interval by Gender",
             x = "Gender",
             y = "Value") +
        theme_minimal() +
        theme(
            plot.title = element_text(size = 15, face = "bold", hjust = 0.5),
            axis.title.x = element_text(size = 12),
            axis.title.y = element_text(size = 12),
            axis.text = element_text(size = 10),
            legend.position = "none"
        ) +
        geom_text(aes(label = round(Mean, 2)), vjust = -0.3, hjust = 1.5, color = "gray30", size = 4) +  # Adicionando valor da média
        geom_text(aes(y = CI_Lower, label = round(CI_Lower, 2)), vjust = 1.5, color = "gray30", size = 4) +  # Valor do limite inferior
        geom_text(aes(y = CI_Upper, label = round(CI_Upper, 2)), vjust = -1, color = "gray30", size = 4) +  # Valor do limite superior
        geom_text(aes(y = SD, label = round(SD, 2)), vjust = -0.3, hjust = -0.5, color = "blue", size = 4) +  # Valor do desvio padrão positivo
        geom_point(aes(y = Mean), size = 4, color = "gray30") +  # Adicionando ponto da média
        geom_point(aes(y = SD), size = 3, color = "blue", shape = 15) +  # Ponto para desvio padrão positivo
        ylim(y_min - abs(0.1 * y_min), y_max + abs(0.1 * y_max))  # Ajustando os limites do eixo y
    
    return(dpr)
}

# Função para calcular a média (ou qualquer outra estatística de interesse)
mean_stat <- function(data, indices) {
    return(mean(data[indices]))
}

# Função para realizar bootstrap e plotar gráfico
bootstrap_ci_plot <- function(data, d_art_col, sex_col, n_bootstrap = 1000, conf_level = 0.95) {
    # Dividir dados por sexo
    data_split <- split(data, data[[sex_col]])
    
    # Lista para armazenar plots
    plots <- list()
    
    for (sex_value in names(data_split)) {
        # Realizar o bootstrap
        boot_results <- boot(data_split[[sex_value]][[d_art_col]], statistic = mean_stat, R = n_bootstrap)
        
        # Calcular intervalo de confiança
        ci <- boot.ci(boot_results, type = "perc", conf = conf_level)
        
        # Extrair estatísticas
        boot_means <- boot_results$t
        original_mean <- mean(data_split[[sex_value]][[d_art_col]])
        lower_ci <- ci$percent[4]
        upper_ci <- ci$percent[5]
        
        # Criar dataframe para ggplot
        df <- data.frame(boot_means)
        
        # Plotar histograma com intervalo de confiança
        plot <- ggplot(df, aes(x = boot_means)) +
            geom_histogram(binwidth = 0.1, fill = "skyblue", color = "black", alpha = 0.7) +
            geom_vline(aes(xintercept = original_mean), color = "red", linetype = "dashed", size = 1, show.legend = TRUE) +
            geom_vline(aes(xintercept = lower_ci), color = "blue", linetype = "dotted", size = 1) +
            geom_vline(aes(xintercept = upper_ci), color = "blue", linetype = "dotted", size = 1) +
            labs(
                title = paste("Bootstrap Distribution with Confidence Interval for", sex_value),
                x = "Bootstrap Means",
                y = "Frequency"
            ) +
            theme_minimal()
        
        # Adicionar plot à lista
        plots[[sex_value]] <- plot
    }
    
    # Organizar os plots um abaixo do outro
    grid_plot <- do.call(grid.arrange, c(plots, ncol = 1))
    
    return(grid_plot)
}

Estatísticas gerais

Para a amostra inteira, com nível de confiança de 95%, temos que a média populacional pode estar descrita no intervalo de confiança abaixo:

print(plot_with_statistics(iat))

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

Estatísticas por gênero

Agrupando por gênero, temos as seguintes observações:

print(plot_with_statistics_by_gender(iat, "d_art", "sex"))

Podemos observar que, em média, as mulheres que participaram do experimento tiveram uma associação implícita (medida pelo IAT) contra a matemárica positiva e moderada (média 0.47, desv. padrão 0.51, N = 77). Homens tiveram uma associação positiva e fraca contra a matemática, portanto menor que a das mulheres (média 0.07, desv. padrão 0.57, N = 24).

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) + 
    geom_text(aes(y = conf.low, label = round(conf.low, 2)), vjust = 1.5, color = "gray30", size = 4) +  # Valor do limite inferior
    geom_text(aes(y = conf.high, label = round(conf.high, 2)), vjust = -0.5, color = "gray30", size = 4) +  # Valor do limite superior
    ylim(-1, 1) + 
    labs(x = "Diferença", 
         y = "IAT mulheres - homens")

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

Houve portanto uma considerável diferença entre homens e mulheres (diferença das médias 0.41, 95% CI [0.16, 0.66]).

A partir desta amostra, estimamos que:

* mulheres têm uma associação positiva consideravelmente mais forte, com uma média que provavelmente está entre 0.35 e  0.58 ponto na escala IAT, o suficiente para diferenciar uma associação neutra de uma considerável contra a matemática.
* homens têm uma associação positiva menos 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.

Estatísticas com bootstrap

Utilizando bootstrap nas visualizações, podemos observar o seguinte:

bootstrap_ci_plot(iat, d_art_col = "d_art", sex_col = "sex")
## 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.

## TableGrob (2 x 1) "arrange": 2 grobs
##   z     cells    name           grob
## m 1 (1-1,1-1) arrange gtable[layout]
## f 2 (2-2,1-1) arrange gtable[layout]

O método escolhido evidencia uma estrutura robusta e flexível para realizar bootstrap em uma ampla variedade de estatísticas e conjuntos de dados. Em comparação com pacotes alternativos, a função boot se destaca pela sua abrangência, integração com outras funções do R e comunidade de usuários ativa.