Aula 4 — Comparação de Médias: Teste t, Testes Não Paramétricos e ANOVA

Objetivos da aula

Nesta aula, vamos organizar uma sequência de análises estatísticas voltadas à comparação de médias em situações comuns em experimentos biológicos e agronômicos. O foco será em três contextos:

comparação entre dois grupos independentes;
comparação entre duas condições medidas nos mesmos indivíduos;
comparação entre três ou mais grupos.

Ao longo da aula, usaremos dados reais, discutindo a lógica estatística, as hipóteses, as fórmulas e a implementação em R.

Pacotes utilizados

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.0     ✔ readr     2.2.0
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.2     ✔ tibble    3.3.1
✔ lubridate 1.9.5     ✔ tidyr     1.3.2
✔ purrr     1.2.1     
── 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(gsheet)
library(epifitter)
library(emmeans)

Warning: pacote 'emmeans' foi compilado no R versão 4.5.3

Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'

library(multcomp)

Warning: pacote 'multcomp' foi compilado no R versão 4.5.3

Carregando pacotes exigidos: mvtnorm

Warning: pacote 'mvtnorm' foi compilado no R versão 4.5.3

Carregando pacotes exigidos: survival
Carregando pacotes exigidos: TH.data

Warning: pacote 'TH.data' foi compilado no R versão 4.5.3

Carregando pacotes exigidos: MASS

Anexando pacote: 'MASS'

O seguinte objeto é mascarado por 'package:dplyr':

    select


Anexando pacote: 'TH.data'

O seguinte objeto é mascarado por 'package:MASS':

    geyser

library(multcompView)

Warning: pacote 'multcompView' foi compilado no R versão 4.5.3

library(rstatix)


Anexando pacote: 'rstatix'

O seguinte objeto é mascarado por 'package:MASS':

    select

O seguinte objeto é mascarado por 'package:stats':

    filter

Parte 1 — Leitura dos dados e cálculo da AUDPC

Leitura da planilha `curve`

O primeiro conjunto de dados contém avaliações de severidade ao longo do tempo, em dois tratamentos de irrigação. Antes de comparar tratamentos, é útil visualizar a progressão temporal da doença. Utilizaremos o pacote gsheet para ler os dados diretamente da nuvem, garantindo a reprodutibilidade sem depender de planilhas locais.

# Carregando os dados via gsheet
dados <- gsheet2tbl("[https://docs.google.com/spreadsheets/d/1bq2N19DcZdtax2fQW9OHSGMR0X2__Z9T/edit?gid=1807247585#gid=1807247585](https://docs.google.com/spreadsheets/d/1bq2N19DcZdtax2fQW9OHSGMR0X2__Z9T/edit?gid=1807247585#gid=1807247585)")

glimpse(dados)

Rows: 60
Columns: 4
$ Irrigation <chr> "Furrow", "Furrow", "Furrow", "Furrow", "Furrow", "Furrow",…
$ rep        <dbl> 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2,…
$ day        <dbl> 0, 0, 0, 7, 7, 7, 14, 14, 14, 21, 21, 21, 28, 28, 28, 35, 3…
$ severity   <dbl> 0.01, 0.01, 0.01, 0.04, 0.04, 0.04, 0.10, 0.10, 0.11, 0.11,…

head(dados)

Resumo gráfico da severidade ao longo do tempo

A seguir, calculamos a média e o desvio-padrão da severidade, em porcentagem, para cada combinação de tratamento e tempo de avaliação.

dados |>
  group_by(Irrigation, day) |>
  summarise(
    mean_sev = mean(severity * 100),
    sd_sev = sd(severity * 100),
    .groups = "drop"
  ) |>
  ggplot(aes(x = day, y = mean_sev, color = Irrigation, group = Irrigation)) +
  geom_point(size = 2) +
  geom_line() +
  geom_errorbar(aes(ymin = mean_sev - sd_sev,
                    ymax = mean_sev + sd_sev),
                width = 0.1) +
  ylim(0, 100) +
  labs(
    x = "Dias após a avaliação inicial",
    y = "Severidade média (%)",
    color = "Irrigação"
  ) +
  theme_minimal()

Área abaixo da curva de progresso da doença

Em estudos epidemiológicos, muitas vezes desejamos resumir toda a curva de progresso em uma única medida. Uma das medidas mais usadas é a AUDPC (Area Under the Disease Progress Curve), ou área abaixo da curva de progresso da doença.

A fórmula geral da AUDPC, pelo método dos trapézios, é:

\[ AUDPC = \sum_{i=1}^{n-1} \left( \frac{y_i + y_{i+1}}{2} \right)(t_{i+1} - t_i) \]

em que: - \(y_i\) e \(y_{i+1}\) são as severidades observadas em duas avaliações consecutivas; - \(t_i\) e \(t_{i+1}\) são os tempos correspondentes; - \(n\) é o número de avaliações.

Essa estatística resume a intensidade da doença ao longo do tempo. Quanto maior a AUDPC, maior a quantidade acumulada de doença no período avaliado.

dados2 <- dados |>
  group_by(Irrigation, rep) |>
  summarise(
    audpc = AUDPC(day, severity),
    .groups = "drop"
  ) |>
  mutate(Irrigation = as.factor(Irrigation))

dados2

Estatísticas descritivas da AUDPC

Antes de qualquer teste formal, é recomendável resumir os dados por tratamento.

resumo_audpc <- dados2 |>
  group_by(Irrigation) |>
  summarise(
    mean_audpc = mean(audpc),
    sd_audpc = sd(audpc),
    n = n(),
    se_audpc = sd_audpc / sqrt(n),
    ic_95 = qt(0.975, df = n - 1) * se_audpc,
    .groups = "drop"
  )

resumo_audpc

Gráfico com intervalo de confiança de 95%

O intervalo de confiança para a média pode ser escrito como:

\[ \bar{x} \pm t_{\alpha/2,\, n-1} \times \frac{s}{\sqrt{n}} \]

em que \(t_{\alpha/2,\, n-1}\) é o quantil da distribuição t de Student com \(n-1\) graus de liberdade. O erro padrão \(s/\sqrt{n}\) reflete a incerteza da média amostral.

resumo_audpc |>
  ggplot(aes(x = Irrigation, y = mean_audpc)) +
  geom_col(width = 0.5, fill = "skyblue", color = "black") +
  geom_point(size = 2) +
  geom_errorbar(aes(ymin = mean_audpc - ic_95,
                    ymax = mean_audpc + ic_95),
                width = 0.05) +
  ylim(0, 20) +
  labs(x = "Irrigação", y = "AUDPC média (IC 95%)") +
  theme_minimal()

Devido à sobreposição dos intervalos de confiança, já podemos inferir visualmente que não há diferença estatística evidente entre os tratamentos.

Parte 2 — Teste t para dois grupos independentes

Quando usar?

O teste t para amostras independentes é usado quando queremos comparar as médias de dois grupos distintos, assumindo independência entre as observações. Neste caso, vamos comparar a AUDPC entre os dois tratamentos de irrigação.

Hipóteses

As hipóteses do teste são:

\[H_0: \mu_1 = \mu_2\] \[H_1: \mu_1 \neq \mu_2\]

ou seja, a hipótese nula afirma que as médias populacionais são iguais, enquanto a hipótese alternativa afirma que elas diferem.

Aplicação em R

Quando as variâncias são assumidas como diferentes, usa-se a versão de Welch do teste t:

teste_t_audpc <- t.test(audpc ~ Irrigation, data = dados2)
teste_t_audpc


    Welch Two Sample t-test

data:  audpc by Irrigation
t = -1.1983, df = 3.1984, p-value = 0.312
alternative hypothesis: true difference in means between group Drip and group Furrow is not equal to 0
95 percent confidence interval:
 -1.4140759  0.6207425
sample estimates:
  mean in group Drip mean in group Furrow 
            13.39333             13.79000

Interpretação

O p-valor deu 0.312. Como o valor de p é maior que o nível de significância de 5% (0.05), não rejeitamos a \(H_0\), ou seja, os tratamentos não diferem estatisticamente, confirmando o que visualizamos no gráfico com os intervalos de confiança.

Parte 3 — Teste t pareado

Situação experimental

Agora usaremos a planilha escala, em que os mesmos avaliadores foram medidos em duas condições de avaliação. Nesse caso, as observações não são independentes, pois cada avaliador aparece nas duas condições. Logo, o teste adequado é o teste t pareado.

escala <- gsheet2tbl("[https://docs.google.com/spreadsheets/d/1bq2N19DcZdtax2fQW9OHSGMR0X2__Z9T/edit?gid=1729131173#gid=1729131173](https://docs.google.com/spreadsheets/d/1bq2N19DcZdtax2fQW9OHSGMR0X2__Z9T/edit?gid=1729131173#gid=1729131173)")

# Visualização inicial
escala |>
  ggplot(aes(x = assessment, y = acuracia)) +
  geom_boxplot(outlier.colour = NA) +
  geom_jitter(width = 0.1, alpha = 0.7) +
  ylim(0, 1) +
  labs(x = "Tipo de avaliação", y = "Acurácia") +
  theme_minimal()

Organização dos dados em formato largo

Para o teste pareado, é conveniente reorganizar a tabela de forma que cada linha corresponda a um avaliador e cada coluna a uma condição.

escala_wide <- escala |>
  dplyr::select(rater, assessment, acuracia) |>
  pivot_wider(
    names_from = assessment,
    values_from = acuracia
  )

escala_wide

Parte 4 — Verificação da normalidade e alternativa não paramétrica

Inspeção gráfica e teste de Shapiro-Wilk

O teste t pareado assume que as diferenças entre pares seguem uma distribuição normal. Podemos checar as condições isoladamente ou direto nas diferenças.

# Checando o grupo Unaided
qqnorm(escala_wide$Unaided)
qqline(escala_wide$Unaided)

shapiro.test(escala_wide$Unaided) # p-value < 0.05 rejeita normalidade


    Shapiro-Wilk normality test

data:  escala_wide$Unaided
W = 0.7748, p-value = 0.007155

# Checando o grupo Aided1
qqnorm(escala_wide$Aided1)
qqline(escala_wide$Aided1)

shapiro.test(escala_wide$Aided1)


    Shapiro-Wilk normality test

data:  escala_wide$Aided1
W = 0.92852, p-value = 0.4335

Como o primeiro conjunto de dados é não-normal e o segundo é normal, devemos seguir com uma análise estatística não-paramétrica.

Teste não paramétrico de Wilcoxon pareado

Se a suposição de normalidade for claramente violada, a alternativa clássica é o teste de postos sinalizados de Wilcoxon. É o teste não-paramétrico equivalente ao t pareado, para variáveis dependentes. (Se fossem independentes, usaríamos o Mann-Whitney).

teste_wilcox <- wilcox.test(
  escala_wide$Unaided,
  escala_wide$Aided1,
  paired = TRUE
)

Warning in wilcox.test.default(escala_wide$Unaided, escala_wide$Aided1, : não é
possível computar o valor de p exato com o de desempate

teste_wilcox


    Wilcoxon signed rank test with continuity correction

data:  escala_wide$Unaided and escala_wide$Aided1
V = 0, p-value = 0.005889
alternative hypothesis: true location shift is not equal to 0

Aplicação do teste t pareado (para fins pedagógicos)

Apenas para ilustrar como rodar o teste t pareado no R:

teste_t_pareado <- t.test(
  escala_wide$Aided1,
  escala_wide$Unaided,
  paired = TRUE
)

teste_t_pareado


    Paired t-test

data:  escala_wide$Aided1 and escala_wide$Unaided
t = 4.4214, df = 9, p-value = 0.001668
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 0.1147647 0.3552353
sample estimates:
mean difference 
          0.235

Visualização das diferenças

No teste pareado, o interesse recai sobre a diferença entre as duas medidas do mesmo indivíduo. Visualizar explicitamente as diferenças entre pares é excelente para compreender a análise:

diferencas <- escala_wide |>
  mutate(diff = Aided1 - Unaided)

diferencas |>
  ggplot(aes(x = diff)) +
  geom_histogram(bins = 10, color = "white", fill = "steelblue") +
  labs(x = "Diferença pareada (Aided1 - Unaided)", y = "Frequência") +
  theme_minimal()

Parte 5 — ANOVA para comparação de três ou mais grupos

Situação experimental

Agora analisaremos a planilha micelial, que contém dados de taxa de crescimento micelial (tcm) para diferentes espécies. Quando o objetivo é comparar a média de uma variável resposta entre três ou mais grupos, a abordagem clássica é a análise de variância (ANOVA).

micelial <- gsheet2tbl("[https://docs.google.com/spreadsheets/d/1bq2N19DcZdtax2fQW9OHSGMR0X2__Z9T/edit?gid=959387827#gid=959387827](https://docs.google.com/spreadsheets/d/1bq2N19DcZdtax2fQW9OHSGMR0X2__Z9T/edit?gid=959387827#gid=959387827)")

set.seed(1) # Padroniza a distribuição dos pontos de jitter
micelial |>
  ggplot(aes(x = reorder(especie, tcm), y = tcm)) + # Ordena do menor pro maior
  geom_boxplot(outlier.colour = NA) +
  geom_jitter(width = 0.2, alpha = 0.7) +
  labs(x = "Espécie", y = "Taxa de crescimento micelial") +
  theme_minimal()

Ajuste do modelo de ANOVA

A hipótese nula (\(H_0\)) da ANOVA é que todas as médias são iguais. A alternativa (\(H_1\)) é que pelo menos uma difere das demais.

m1 <- lm(tcm ~ especie, data = micelial)
anova(m1)

Como o p-valor é muito pequeno (\(p < 0.05\)), rejeitamos \(H_0\). Pelo menos um grupo difere dos demais.

Parte 6 — Diagnóstico dos pressupostos da ANOVA

A ANOVA clássica assume normalidade dos resíduos e homogeneidade de variâncias (homocedasticidade).

# 1. Normalidade dos Resíduos
hist(m1$residuals)

qqnorm(m1$residuals)
qqline(m1$residuals)

shapiro.test(m1$residuals)


    Shapiro-Wilk normality test

data:  m1$residuals
W = 0.9821, p-value = 0.8782

# p-valor grande, não rejeita H0, os resíduos têm distribuição normal.

# 2. Homogeneidade de variâncias (Teste de Bartlett)
bartlett.test(tcm ~ especie, data = micelial)


    Bartlett test of homogeneity of variances

data:  tcm by especie
Bartlett's K-squared = 4.4367, df = 4, p-value = 0.3501

# p-valor grande, não rejeita H0, as variâncias são homogêneas.

Parte 7 — Médias ajustadas com `emmeans` e comparação de Tukey

Depois que a ANOVA indica diferença global entre médias, surge a pergunta: quais grupos diferem entre si? Fazer vários testes t independentes inflaciona a probabilidade de erro do tipo I. O método correto é usar testes de comparações múltiplas como o Tukey.

# Médias marginais estimadas (ajustadas pelo modelo)
medias_m1 <- emmeans(m1, ~ especie)
medias_m1

 especie emmean     SE df lower.CL upper.CL
 Fasi     1.572 0.0559 25    1.457     1.69
 Faus     1.237 0.0559 25    1.122     1.35
 Fcor     1.322 0.0559 25    1.207     1.44
 Fgra     0.912 0.0559 25    0.797     1.03
 Fmer     1.427 0.0559 25    1.312     1.54

Confidence level used: 0.95

# Comparações múltiplas de Tukey
pairs(medias_m1, adjust = "tukey")

 contrast    estimate    SE df t.ratio p.value
 Fasi - Faus    0.335 0.079 25   4.241  0.0023
 Fasi - Fcor    0.250 0.079 25   3.165  0.0302
 Fasi - Fgra    0.660 0.079 25   8.356 <0.0001
 Fasi - Fmer    0.145 0.079 25   1.836  0.3765
 Faus - Fcor   -0.085 0.079 25  -1.076  0.8169
 Faus - Fgra    0.325 0.079 25   4.115  0.0031
 Faus - Fmer   -0.190 0.079 25  -2.405  0.1469
 Fcor - Fgra    0.410 0.079 25   5.191  0.0002
 Fcor - Fmer   -0.105 0.079 25  -1.329  0.6761
 Fgra - Fmer   -0.515 0.079 25  -6.520 <0.0001

P value adjustment: tukey method for comparing a family of 5 estimates

Letras de agrupamento

Uma forma prática de resumir os resultados é usar letras. Médias seguidas pela mesma letra não diferem entre si ao nível de significância adotado.

medias_m1b <- cld(medias_m1, Letters = letters)
medias_m1b

Gráfico das médias com intervalos de confiança e letras

medias_m1b |>
  ggplot(aes(x = reorder(especie, emmean), y = emmean, label = .group)) +
  geom_point(size = 2) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.06) +
  geom_text(hjust = 0, nudge_y = 0.12, size = 4) +
  ylim(0.5, 2) +
  coord_flip() +
  labs(
    x = "",
    y = "Taxa de crescimento micelial (mm/dia)"
  ) +
  theme_classic()

Parte 8 — Síntese conceitual da aula

1. Teste t para amostras independentes Usado para comparar a média de dois grupos independentes. - Unidade analítica em cada grupo é diferente. - Exemplo desta aula: comparação da AUDPC entre tratamentos de irrigação.

2. Teste t pareado Usado quando as duas medidas vêm dos mesmos indivíduos ou de unidades emparelhadas. - O teste é realizado sobre as diferenças entre pares. - Exemplo desta aula: comparação da acurácia dos mesmos avaliadores sob duas condições.

3. Teste de Wilcoxon pareado Alternativa não paramétrica ao teste t pareado. - Útil quando a distribuição das diferenças é muito assimétrica ou quando há clara violação da normalidade.

4. ANOVA Usada para comparar três ou mais médias simultaneamente. - Primeiro testa-se o efeito global. - Depois, se necessário, investigam-se comparações múltiplas.

5. emmeans e Tukey Permitem estimar médias ajustadas e comparar grupos com controle do erro associado a múltiplas comparações.

Parte 9 — Sugestões de interpretação para o relatório

Ao redigir resultados, procure separar claramente:

Descrição dos dados, com médias e medidas de dispersão;
Resultado do teste global, com estatística e valor de p;
Comparações específicas, quando cabíveis;
Interpretação biológica dos resultados.

Exemplo de estrutura textual: A AUDPC média foi maior em um dos tratamentos de irrigação. O teste t para amostras independentes indicou diferença significativa entre as médias dos grupos (p < 0,05). No experimento com avaliadores, a acurácia diferiu entre as condições quando analisada de forma pareada. Para os dados de crescimento micelial, a ANOVA indicou efeito de espécie, e as comparações de Tukey permitiram identificar quais espécies apresentaram médias distintas.

Objetivos da aula

Pacotes utilizados

Parte 1 — Leitura dos dados e cálculo da AUDPC

Leitura da planilha curve

Resumo gráfico da severidade ao longo do tempo

Área abaixo da curva de progresso da doença

Estatísticas descritivas da AUDPC

Gráfico com intervalo de confiança de 95%

Parte 2 — Teste t para dois grupos independentes

Quando usar?

Hipóteses

Aplicação em R

Interpretação

Parte 3 — Teste t pareado

Situação experimental

Organização dos dados em formato largo

Parte 4 — Verificação da normalidade e alternativa não paramétrica

Inspeção gráfica e teste de Shapiro-Wilk

Teste não paramétrico de Wilcoxon pareado

Aplicação do teste t pareado (para fins pedagógicos)

Visualização das diferenças

Parte 5 — ANOVA para comparação de três ou mais grupos

Situação experimental

Ajuste do modelo de ANOVA

Parte 6 — Diagnóstico dos pressupostos da ANOVA

Parte 7 — Médias ajustadas com emmeans e comparação de Tukey

Letras de agrupamento

Gráfico das médias com intervalos de confiança e letras

Parte 8 — Síntese conceitual da aula

Parte 9 — Sugestões de interpretação para o relatório

Leitura da planilha `curve`

Parte 7 — Médias ajustadas com `emmeans` e comparação de Tukey