A questão

O experimento

Obter os dados experimentais

library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.2.1
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.8     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.1
## ✔ readr   2.1.2     ✔ forcats 0.5.2
## Warning: package 'tibble' was built under R version 4.2.1
## Warning: package 'stringr' was built under R version 4.2.1
## Warning: package 'forcats' was built under R version 4.2.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
mosquitos <- read_csv2("mosquitos.csv")
## ℹ Using "','" as decimal and "'.'" as grouping mark. Use `read_delim()` for more control.
## Rows: 25 Columns: 2── Column specification ────────────────────────────────────────────────────────
## Delimiter: ";"
## dbl (2): agua, cerveja
## ℹ 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.
summary(mosquitos)
##       agua          cerveja    
##  Min.   :12.00   Min.   :17.0  
##  1st Qu.:16.50   1st Qu.:20.0  
##  Median :20.00   Median :24.0  
##  Mean   :19.22   Mean   :23.6  
##  3rd Qu.:22.00   3rd Qu.:27.0  
##  Max.   :24.00   Max.   :31.0  
##  NA's   :7
mosquitos %>% gather('agua', 'cerveja', key = "bebida", value = "mordidas", na.rm = TRUE) -> zica
table(zica)
##          mordidas
## bebida    12 13 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 31
##   agua     1  1  2  1  0  1  2  2  2  3  1  2  0  0  0  0  0  0
##   cerveja  0  0  0  0  1  1  2  4  3  0  1  3  1  1  3  2  1  2

teste-t (duas maneiras)

tt1 <- t.test(mordidas ~ bebida, data = zica, var.equal = FALSE)
tt1
## 
##  Welch Two Sample t-test
## 
## data:  mordidas by bebida
## t = -3.6582, df = 39.113, p-value = 0.0007474
## alternative hypothesis: true difference in means between group agua and group cerveja is not equal to 0
## 95 percent confidence interval:
##  -6.798084 -1.957472
## sample estimates:
##    mean in group agua mean in group cerveja 
##              19.22222              23.60000
tt2 <- t.test(mosquitos$cerveja, mosquitos$agua, 
              alternative = "greater", 
              var.equal = FALSE, 
              na.action = "na.omit")
tt2
## 
##  Welch Two Sample t-test
## 
## data:  mosquitos$cerveja and mosquitos$agua
## t = 3.6582, df = 39.113, p-value = 0.0003737
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  2.36165     Inf
## sample estimates:
## mean of x mean of y 
##  23.60000  19.22222
(ggplot(zica, aes(bebida, mordidas))
    + geom_boxplot()
    + stat_sum(colour ="darkgray", alpha = 0.5)
    + scale_size(breaks = 1:2, range = c(3,6))
)

permutação

set.seed(999) ## para reproducibilidade
nsim <- 10000
res <- numeric(nsim) ## reserva espaço para os resultados das diferenças
for (i in 1:nsim) {
    ## padrão: "mistura" as respostas
    perm <- sample(nrow(zica))
    bdat <- transform(zica, mordidas = mordidas[perm])
    ## calcula e guarda as diferenças entre as médias
    res[i] <- mean(bdat[bdat$bebida == "agua", "mordidas"]) -
              mean(bdat[bdat$bebida == "cerveja", "mordidas"])
}

obs <- zica %>% group_by(bebida) %>% 
                summarize(mordidas = mean(mordidas)) %>% 
                pull(mordidas) %>%  ## extrai uma coluna
                diff()          ## diferença entre elementos

Uma outra maneira: obs <- colMeans(mosquitos, na.rm = TRUE) %>% diff()

Anexa o valor observado à lista de resultados das permutações

res <- c(res, obs)

“p-valor” do teste. A proporção dos casos onde res >= obs.

mean(res >= obs)
## [1] 0.00079992

Histograma

hist(res, 
     col = "gray", 
     las = 1,
     main = "permutação 10.000 vezes", 
     xlab = "diferenças", 
     ylab = "Contagem")
abline(v = obs, 
       col = "red", 
       lwd = 3)

Outro tipo de gráfico

par(las = 1,bty = "l")
plot(prop.table(table(round(res, 2))),
     main = "Permutação 50.000 vezes", 
     xlab = "diferenças", 
     ylab = "Proporção", 
     axes = FALSE)
axis(side = 2)
points(obs, 0, pch = 16, cex = 2.0, col = "red")

Teste de randomização de Fisher para duas amostras (médias)

Testes de randomização.

Em 1935, R.A. Fisher introduziu a ideia do teste de randomização, que se baseia em tentar responder à pergunta:
“Será que o padrão observado aconteceu por acaso, ou o padrão indica que a hipótese nula não é verdadeira?”
Um teste de randomização funciona simplesmente enumerando todos os possíveis resultados sob a hipótese nula e, em seguida, registrando onde o resultado observado se encaixa.

Um teste de randomização também é chamado de teste de permutação, pois envolve permutar as observações durante o procedimento de enumeração.

Um p-valor superior unilateral é calculado como a proporção de vezes que as diferenças das médias (ou medianas) na distribuição de permutação são maiores ou iguais à diferença observada.

library(EnvStats)
## Warning: package 'EnvStats' was built under R version 4.2.1
## 
## Attaching package: 'EnvStats'
## The following objects are masked from 'package:stats':
## 
##     predict, predict.lm
## The following object is masked from 'package:base':
## 
##     print.default
tt3 <- twoSamplePermutationTestLocation(mosquitos$cerveja, na.omit(mosquitos$agua), 
                                        fcn = "mean", 
                                        alternative = "greater", 
                                        n.permutations = 10000,
                                        seed = 999)
plot(tt3)

Agora, um teste utilizando Bootstrap

Observa-se duas amostras independentes x e y.
Deseja-se testar a hipótese básica de médias iguais.

\(x = \langle x_{1}, x_{2}, \dots, x_{n}\rangle \rightarrow \bar{x}\)
\(y = \langle y_{1}, y_{2}, \dots, y_{m}\rangle \rightarrow \bar{y}\)
\(dif = \bar{x} - \bar{y}\)

Teste de hipótese:
- \(H_0\): \(\mu_x = \mu_y\) (hipótese nula, ou básica)
- \(H_1\): \(\mu_x > \mu_y\) (hipótese alternativa).

Hipótese e nível alfa:
• \(H_0\): Ambas as amostras da mesma população
• \(H_1\): Ambas as amostras não são da mesma população e \(\mu_x > \mu_y\)
• \(\alpha\) = 0,05

O procedimento de Bootstrap:
- Passo 0: Juntar as duas amostras observadas* em um único conjunto com (n + m) observações
- Passo 1: Obter uma amostra bootstrap de (n + m) observações com reposição do conjunto total
- Passo 2: Calcular a média das n primeiras novas observações e chamar de \(\bar{x}\)
Calcular a média das m novas observações restantes e chamar de \(\bar{y}\)a
Calcular a estatísitica de teste \(t^* = \bar{x} - \bar{y}\)

  • Passo 3: Repetir os passos 1 e 2 um número B (e.g., 10.000) de vezes, obtendo B valores da estatística de teste \(t^*\)
  • Passo 4: O p-valor desejado é estimado como:
    p−valor \(\approx\) (número de vezes {\(t^*\) > dif}) / B

Rejeitar \(H_0\) se p−valor < \(\alpha\)

Notas
• Não se assume a normalidade das populações
• Portanto, a distribuição das amostras (i.e., a distribuição de todas as possíveis estatísiticas \(t^*\)) não segue, em geral, a distribuição t de Student

em R

library(boot)
cerveja <- mosquitos$cerveja
agua <- na.omit(mosquitos$agua)
set.seed(999)
library(MKinfer)
## Warning: package 'MKinfer' was built under R version 4.2.1
(tt4 <- boot.t.test(cerveja, agua,
                   alternative = 'greater',
                   R = 10000))
## 
##  Bootstrapped Welch Two Sample t-test
## 
## data:  cerveja and agua
## bootstrapped p-value = 6e-04 
## 95 percent bootstrap percentile confidence interval:
##  2.508889      Inf
## 
## Results without bootstrap:
## t = 3.6582, df = 39.113, p-value = 0.0003737
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  2.36165     Inf
## sample estimates:
## mean of x mean of y 
##  23.60000  19.22222

p-valor

tt4$p.value
## [1] 0.000373701