Olá!

Neste post pretendo explicar como realizar um teste de hipóteses usando o R, e ainda também falar a respeito do que é o p-valor que tanto aparece em artigos científicos.

ToothGrowth Data set

Neste data set imbutido na linguagem (lembre-se de carregar o tidyverse) chamado ToothGrowth encontramos 60 observações dividas em três variáveis:

len
supp
dose

Esses dados foram obtidos a partir de uma pesquisa, bem melhor explicada neste post, mas que no geral demonstra o efeito da vitamina C no crescimento de odontoblastos (células responsáveis pelo crescimento dentário) em diversos porquinhos-da-índia, com diferentes doses.

ggplot(ToothGrowth) +
  geom_freqpoly(aes(x = len, color = dose), binwidth = 1) +
  theme_classic()

ggplot(ToothGrowth) + geom_freqpoly(aes(x = len, color = dose), binwidth = 1) + theme_classic()

Parece que existe uma clara influência da dose de vitamina aplicada nos porquinhos-da-índia no crescimento de seus odontoblastos, porém precisamos ser céticos e obter algumas evidências a respeito do que estamos clamando.

Primeiro observamos as médias

Para começar a compreender a relação da dose de vitamina com o crescimento dentário, podemos observar as médias usando a função summarise() do pacote dplyr.

ToothGrowth %>%
  group_by(dose) %>%
  summarise(media = mean(len))

# A tibble: 3 × 2
  dose  media
  <fct> <dbl>
1 0.5    10.6
2 1      19.7
3 2      26.1

Parece que existe uma diferença entre as médias de cada dose de vitamina aplicada, vamos analisar as diferenças observáveis entre as médias:

diferenca_observavel1 <- 19.735 - 10.605
diferenca_observavel2 <- 26.100 - 19.735
diferenca_observavel3 <- 26.100 - 10.605

Chegando até aqui já respondemos algumas perguntas:

Seriam os dentes dos porquinhos-da-índia que receberam uma dose maior de vitamina C, geralmente maiores?

#uma outra forma de se obter a variação por dose

dose_05 <- ToothGrowth %>%
              filter(dose == 0.5)

dose_1 <- ToothGrowth %>%
              filter(dose == 1)
              
dose_2 <- ToothGrowth %>%
              filter(dose == 2)

Podemos checar se isso é verdade da seguinte forma:

#Dentro de diferenca_observavel1 temos a diferença entre a média das doses 0.5 e 1

mean(dose_1$len) - mean(dose_05$len) == diferenca_observavel1

[1] TRUE

Quais os maiores odontoblastos de cada grupos? E os menores?

ToothGrowth %>%
  group_by(dose) %>%
  summarise(minimo = min(len), maximo = max(len))

# A tibble: 3 × 3
  dose  minimo maximo
  <fct>  <dbl>  <dbl>
1 0.5      4.2   21.5
2 1       13.6   27.3
3 2       18.5   33.9

Então é isso, parece que a dose da vitamina C tem efeito direto sobre o crescimento do dente dos porquinhos-da-índia, na verdade não… o que realmente acontece é que se repetíssemos o experimento com diferentes grupos de porquinhos-da-índia, provavelmente obteríamos uma resposta diferente daquela que observamos previamente.

Por essa razão, essas variáveis são conhecidas como variáveis aleatórias, elas podem conter diversos valores diferentes.

Ou seja, todas as vezes que repetirmos o experimento, obteríamos diferentes médias para os grupos que possuimos (dose), sendo assim a variável len é considerada aleatória.

Vamos entender melhor como funcionam as variáveis aleatórias.

Variáveis aleatórias

Vamos compreender como as médias mudam quando pegamos amostras do nosso data set:

controle <- ToothGrowth$len %>% sample(12)

Podemos fazer isso para construir uma distribuição completamente aleatória a partir dos nossos dados (pois não temos como repetir o experimento, de fato, então modelamos como seria a população caso não houvesse aplicação de dose alguma).

Perceba como a média muda a cada amostragem, repita o comando de amostragem em seu R usando a função sample().

media_controle <- ToothGrowth$len %>% sample(12) %>% mean()

print(media_controle)

media_controle <- ToothGrowth$len %>% sample(12) %>% mean()

print(media_controle)

[1] 19.875

[1] 15.39167

Perceba como a média varia a cada medida, podemos fazer isso várias vezes para aprender ainda mais sobre a distribuição dessa variável aleatória.

Conforme dito anteriormente precisamos nos manter céticos, para isso partimos do pressuposto de que aleatoriamente a diferença observada poderia não significar muito, ou seja, poderia ser algo que aconteceu ao acaso.

Para obter evidências e minimizar essa incerteza, usamos a hipótese nula.

Hipótese Nula

Vamos voltar a analisar as diferenças observáveis de antes. Como precisamos ser céticos temos de nos perguntar:

Como podemos afirmar que as diferenças observáveis acontecem por causa da dose aplicada?
O que aconteceria se dessemos a todos os porquinhos-da-índia a mesma dose?
Veríamos diferenças tão grandes quanto essas?

Neste momento estamos trabalhando com a Hipótese Nula, cujo nome vem para nos lembrar de que estamos agindo como céticos: damos crédito a possibilidade de que pode não existir nenhuma diferença verdadeira em nossas observações, e que elas podem ter acontecido ao acaso.

Nula, pois nos perguntamos o que aconteceria com a diferença observável caso não houvesse nenhum tipo de aplicação acontecendo no experimento, caso nada tivesse acontecido, o que aconteceria?

Construindo um grupo controle para testar a hipótese nula

Podemos observar, quantas vezes quisermos, as diferenças entre os tamanhos quando não estamos comparando as doses. Podemos fazer isso ao amostrar um número aleatório de porquinhos-da-índia, com nenhuma dose aplicada, a partir de uma população aleatória que contenha quaisquer 12 porquinhos-da-índia da população original, e dai gravar o valor da diferença entre a média de dois grupos aleatoriamente. O processo escrito em R, aprendi com o prof. Rafa, no livro “Data analysis for the life sciences”:

media_controle <- ToothGrowth$len %>% sample(12) %>% mean()

media_tratamento <- ToothGrowth$len %>% sample(12) %>% mean()

print(media_tratamento - media_controle)

Inclusive é possível fazer isso diversas vezes usando loops no R:

n <- 10000 
null <- vector("numeric",n) #primeiro criamos um vetor com 10000 espaços
for (i in 1:n) { #looping para criar um vetor nulo
media_controle <- ToothGrowth$len %>% sample(12) %>% mean()
media_tratamento <- ToothGrowth$len %>% sample(12) %>% mean()
null[i] <- mean(media_tratamento) - mean(media_controle)
}

print(null)

Os valores em null são o que chamamos de distribuição nula.

Então, qual a porcentagem de null é maior ou igual do que as diferenças observáveis?

mean(null >= diferenca_observavel1)

mean(null >= diferenca_observavel2)

mean(null >= diferenca_observavel3)

[1] 5e-04

[1] 0.0111

[1] 0

Interpretando os resultados

Parece que apenas uma pequena porcentagem dentre as nossas 10.000 simulações é igual ou maior as diferenças observadas na hipótese nula. Como céticos o que concluimos? Quando não existe efeito de doses, podemos observar uma diferença tão grande quanto observamos antes apenas 0,0001% das vezes (no máximo com a diferenca_observavel1). Esse valor é conhecido como p-value ou p-valor.

P-value ou P-valor

Fizemos uma simulação para contruir o data set null, vamos observar as diferenças observaveis caso não houvessem doses aplicadas:

hist(null, freq=TRUE)

Podemos observar também quais valores são tão maiores quanto as nossas diferenças observadas:

hist(null, freq=TRUE)

abline(v=diferenca_observavel1, col="red", lwd=2)

abline(v=diferenca_observavel2, col="blue", lwd=2)

abline(v=diferenca_observavel3, col="yellow", lwd=2)

Os valores maiores ou iguais as diferenças observáveis são relativamente raros. O que nos garante maior evidência a respeito da nossa hipótese de que a quantidade de vitamina C afeta o crescimento dos odontoblastos.

Foi possível perceber que a maior diferença, entre as doses 0.5 e 2, nem aparecem em nossa distribuição nula, mostrando que essa diferença observável não ocorreu ao acaso e que a hipótese nula pode ser refutada.

Mas ainda assim devemos tomar cuidado!

Grupo controle a partir de outras doses

Para criar um grupo controle com apenas uma dose, podemos filtrar o data set original

ToothGrowth_dot5 <- ToothGrowth %>%
  filter(dose == 0.5)

Agora podemos criar nossos data sets para construir a hipótese nula:

controle_ToothGrowth_dot5 <- ToothGrowth_dot5$len %>% sample(12) %>% mean()

tratamento_ToothGrowth_dot5 <- ToothGrowth_dot5$len %>% sample(12) %>% mean()

print(tratamento_ToothGrowth_dot5 - controle_ToothGrowth_dot5)

[1] -1.483333

Aqui repetimos 10.000 vezes a diferença entre as médias para criar uma distribuição nula:

n <- 10000 
null_dose <- vector("numeric",n) #primeiro criamos um vetor com 10000 espaços
for (i in 1:n) { #looping para criar um vetor nulo
controle_ToothGrowth_dot5 <- ToothGrowth_dot5$len %>% sample(12) %>% mean()
tratamento_ToothGrowth_dot5 <- ToothGrowth_dot5$len %>% sample(12) %>% mean()
null_dose[i] <- mean(tratamento_ToothGrowth_dot5) - mean(controle_ToothGrowth_dot5)
}

print(null)

hist(null_dose, freq=TRUE)

Usando o pacote infer()

O pacote infer nos ajuda a testar hipóteses de uma forma bem fácil e direta Instale e carregue o pacote:

install.packages("infer")

library(infer)

Para calcular a estatística da diferença observada:

F_hat <- ToothGrowth %>% 
  specify(len ~ dose) %>%
  calculate(stat = "F")

Gerando a distribuição nula:

null_dist <- ToothGrowth %>% 
   specify(len ~ dose) %>%
   hypothesize(null = "independence") %>%
   generate(reps = 1000, type = "permute") %>%
   calculate(stat = "F")

Visualizando a diferença observada com a distribuição nula:

visualize(null_dist) +
  shade_p_value(obs_stat = F_hat, direction = "greater")

Falaremos ainda mais sobre isso em outros posts quando usaremos todo esse conhecimento aplicado a um projeto! Até mais!

Teste de hipótese e p-valor