Teste de Kruskal-Wallis para comparação de 3 ou mais populações independentes

EST212 - Bioestatística

Helgem de Souza

Introdução

Antes de iniciar nossa aula:

  1. Abra o RStudio.

  2. Crie um novo script.

  3. Crie uma pasta na área de trabalho com seu nome.

  4. Defina a pasta como diretório de trabalho

  5. Salve o script criado na pasta com o nome "aula16_est212.R".

  6. Baixe do Moodle o arquivo penguins.csv , .

  7. Leia o arquivo listado acima para o objetos penguin.

Comparação de 3 ou mais Populações Independentes

Conforme vimos anteriormente, muitas vezes é necessária a comparação de 3 ou mais populações independentes. Para isso, aprendemos a realizar a Análise de Variância (ANOVA).

A ANOVA tinha alguns pressupostos fundamentais, dos quais destacaremos dois:

  • Os resíduos devem apresentar normalidade, salvo o caso de grandes amostras (grupos com mais de 30 elementos)

  • As variâncias dos grupos devem ser estatisticamente iguais, ou seja, deve ocorrer a homogeneidade de variâncias.

Quando esses requisitos não são cumpridos, não podemos aplicar as conclusões obtidas de forma estatisticamente correta.

Vejamos um exemplo.

Exemplo 1 - Compimento de bicos de penguins por ilha

Vamos utilizar o banco de dados de penguins, que usamos no estudo da ANOVA de um fator. Queremos verificar se os comprimentos de bico diferem entre as ilhas.

Primeiramente, vamos analisar visualmente os dados por meio de boxplots.

Exemplo 1 - Compimento de bicos de penguins por ilha

#Boxplot do comprimento de bico por ilha
boxplot(comprimento_bico ~ ilha, data = penguins, 
        xlab = "Ilha", ylab = "Comprimento do bico", 
        col = c("red", "blue", "green"))

É possível perceber que aparentemente existem diferenças tanto nas medianas quanto na variabilidade. Vamos verificar se podemos utilizar a ANOVA.

Exemplo 1 - Compimento de bicos de penguins por ilha

O primeiro passo é verificar se as variâncias são homogêneas. Utilizaremos o teste de Levene, por meio da função Levenetest, do pacote car. As hipóteses são:

  • \(H_0\): Todos os grupos apresentam variâncias iguais

  • \(H_1\): Existe pelo menos um grupo com variância diferente das demais.

#Carregamento do pacote car
library(car)

#Teste de Levene para homogeneidade de variâncias
leveneTest(comprimento_bico ~ ilha, data = penguins)
Levene's Test for Homogeneity of Variance (center = median)
       Df F value    Pr(>F)    
group   2  22.243 8.384e-10 ***
      339                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

O teste de Levene indica que as variâncias não são iguais, ou seja, um dos requisitos da ANOVA não é atendido.

Exemplo 1 - Compimento de bicos de penguins por ilha

Como existe pelo menos um grupo com variância diferente dos demais, não podemos prosseguir com a ANOVA. Mas vamos seguir com a análise para verificar os resíduos:

#ANOVA - Comprimento de bico por ilha
anova_bico <- aov(comprimento_bico ~ilha, data = penguins)
shapiro.test(anova_bico$residuals)

    Shapiro-Wilk normality test

data:  anova_bico$residuals
W = 0.98889, p-value = 0.01042

Note que os resíduos não são normais, ou seja, os próprios comprimentos de bico por ilha não apresentam distribuição normal.

Neste caso, não podemos utilizar os resultados da ANOVA e é necessária a aplicação de um teste em que esses requisitos não sejam mandatórios.

Teste de Kruskal-Wallis

Quando os requisitos para o uso da ANOVA de um fator não são atingidos, o substituto imediato é o Teste de Kruskal-Wallis.

Seu objetivo é semelhante ao da ANOVA. Enquanto a ANOVA verifica se três ou mais amostras possuem a mesma média, o teste de Kruskal-Wallis verifica se três ou mais amostras possuem a mesma mediana. Suas hipóteses são:

  • \(H_0\): As amostras tem a mesma mediana

  • \(H_1\): Pelo menos uma dass amostras apresenta mediana diferente das demais.

Vamos entender lógica da construção do teste.

Teste de Kruskal-Wallis

Considere as seguintes k amostras, com respectivos tamanhos \(n_1, n_2, \cdots, n_k\):

Amostra 1 \(X_{11}\) \(X_{12}\) \(\cdots\) \(X_{1n_1}\)
Amostra 2 \(X_{21}\) \(X_{22}\) \(\cdots\) \(X_{2n_2}\)
\(\vdots\)
Amostra k \(X_{k1}\) \(X_{k2}\) \(X_{kn_k}\)

Considerando a hipótese nula verdadeira, espera-se que quando ranquearmos todas as observações, as amostras apresentem uma soma de ranks semelhante.

Ou seja, se todas as amostras apresentam rank semelhante, a diferença entre o valor observado e o valor esperado seria pequeno.

A seguir apresentaremos a formulaçã matemática, a título de ilustração. Detalhes podem ser obtidos em KVAM; VIDAKOVIC; KIM (2022), páginas 154 e 155.

Teste de Kruskal-Wallis - Formulação

Seja:

  • \(R_i = \sum\limits_{j = 1}^{n_i}r(X_{ij})\) a soma dos ranks da amostra \(i\).

  • \(n = n_1 + n_2 + \cdots + n_k\) o número de elementos de todas as amostras.

É possível demonstrar que a estatística

\[ H = \dfrac{12}{n(n+1)}\sum\limits_{i = 1}^k\dfrac{1}{n_i}\left[R_i - \dfrac{n_i(n+1)}{2}\right]^2 \]

Segue uma distribuição \(\chi^2\) com k-1 graus de liberdade.

A estatística acima é uma medida entre as somas de ranks observados e esperados em cada amostra. Com base nessa distribuição, calcula-se os p-valores para o teste de Kruskal-Wallis.

Como o teste já está implementado no R, não precisaremos realizar os cálculos.

Teste de Kruskal-Wallis no R

O teste de Kruskal-Wallis está implementado no R, por meio da função kruskal.test(formula, data = dados), em que

  • formula - novamente nossa fórmula habitual, do tipo variável ~ grupo

  • dados - representa o conjunto de dados em estudo.

Sua aplicação no R é bastante semelhante à ANOVA. Entretanto, no caso do teste de Kruskal-Wallis não é necessário utilizar a função summary para acessar o p-valor. Ele é exibido por padrão.

Vamos agora concluir o exemplo dos penguins com o teste de Kruskal-Wallis

Exemplo 1 - Compimento de bicos de penguins por ilha

Nos slides anteriores vimos que os dados não apresentam variância constante e nem normalidade. Por esse motivo a ANOVA não é um teste adequado.

Agora que conhecemos o teste adequado, vamos concluir nossa análise por meio da utilização do teste de Kruskal-Wallis, ao nível de 5% de significância:

#Teste de Kruskal-Wallis para comparação de populações
kruskal.test(comprimento_bico ~ ilha, data = penguins)

    Kruskal-Wallis rank sum test

data:  comprimento_bico by ilha
Kruskal-Wallis chi-squared = 52.438, df = 2, p-value = 4.103e-12

Como o p-valor é menor que o nível de significância, podemos afirmar, ao nível de 5% de significância, que em pelo menos uma das ilhas os penguins apresentam um comprimento de bico maior.

Exercício Prático 1

Vamos verificar, ao nível de 5% de significância, se os comprimentos de bico são diferentes por espécie de penguim.

Exercício Prático 1

Inicialmente, vamos fazer uma análise descritiva dos dados

#Mediana do comprimento de bico por espécie
aggregate(comprimento_bico ~ especie, data = penguins, FUN = median)
    especie comprimento_bico
1    Adelie            38.80
2 Chinstrap            49.55
3    Gentoo            47.30

Os penguins da espécie Adelie presentes na amostra apresentam um comprimento mediano de bico inferior, se comparado aos demais.

Antes de realizar o teste para confirmar a tendência ou não, vamos analisar os boxplots.

Exercício Prático 1

# Boxplot - Comprimento de bico por espécie
boxplot(comprimento_bico ~ especie, data = penguins,
        xlab = "Espécie", ylab = "Comprimento de bico", col = "lightblue")

Podemos verificar que os penguins da espécie Adelie tem um bico mais curto, se comparado às demais raças. Vamos executar o teste para verificar se a diferença é significativa.

Exercício Prático 1

Assim como no exemplo, vamos usar a função kruskal.test. As hipóteses são as seguintes:

  • \(H_0:\) Todas as espécies apresentam o mesmo comprimento mediano de bico.

  • \(H_1:\) Pelo menos uma das espécies apresenta comprimento mediano de bico diferente das demais.

#Teste de Kruskal-Wallis - Comprimento de bico por espécie
kruskal.test(comprimento_bico ~ especie, data = penguins)

    Kruskal-Wallis rank sum test

data:  comprimento_bico by especie
Kruskal-Wallis chi-squared = 244.14, df = 2, p-value < 2.2e-16

Como o p-valor é menor que o nível de significância, temos evidências que nos levam a não aceitar a hipótese nula. Ou seja, ao nível de 5% de significância, podemos afirmar que pelo menos uma das espécies apresenta comprimento de bico diferente das demais.

Comparações múltiplas

Apesar da evidência visual, somente podemos confirmar as diferenças por meio de testes de comparações múltiplas.

Trataremos desse assunto na próxima aula.

Referências

KVAM, P.; VIDAKOVIC, B.; KIM, S. Nonparametric statistics with applications to science and engineering with r. [s.l.] John Wiley & Sons, 2022.