Aula 5 - Análise de Dados de Contagem: Da ANOVA ao GLM
1. Introdução e Análise Exploratória
Nesta análise, vamos investigar a eficácia de diferentes inseticidas (spray) na contagem de insetos sobreviventes (count). Este é um exemplo clássico de dados de contagem (números inteiros não negativos).
Primeiro, carregamos os pacotes necessários e realizamos uma análise visual exploratória para entender a distribuição e a variabilidade dos dados.
# Carregar dados nativos do Rinsects <- InsectSprays# Carregar pacoteslibrary(ggplot2)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
✔ lubridate 1.9.5 ✔ tibble 3.3.1
✔ purrr 1.2.1 ✔ tidyr 1.3.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
# Análise visual exploratória# Combinamos jitter (para ver cada ponto individualmente sem sobreposição) # e boxplot (para visualizar medianas e quartis).insects |>ggplot(aes(x = spray, y = count)) +geom_jitter(width =0.1, alpha =0.5) +geom_boxplot(outlier.colour =NA, fill =NA, color ="blue") +theme_classic() +labs(title ="Contagem de insetos por tipo de Spray",x ="Tipo de Spray",y ="Contagem")
2. Modelo Linear Tradicional (ANOVA) e Pressupostos
A abordagem mais comum (e muitas vezes inicial) é ajustar um modelo linear clássico (ANOVA). O modelo linear (\(Y = \beta_0 + \beta_1X + \epsilon\)) assume que os resíduos (\(\epsilon\)) têm: 1. Normalidade: Seguem uma distribuição Gaussiana. 2. Homocedasticidade: Possuem variância constante entre os tratamentos.
Vamos ajustar o modelo e testar essas suposições matematicamente e visualmente.
# Ajuste do modelo linearm1 <-lm(count ~ spray, data = insects)# Resumo do modelo e quadro da ANOVAsummary(m1)
Call:
lm(formula = count ~ spray, data = insects)
Residuals:
Min 1Q Median 3Q Max
-8.333 -1.958 -0.500 1.667 9.333
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.5000 1.1322 12.807 < 2e-16 ***
sprayB 0.8333 1.6011 0.520 0.604
sprayC -12.4167 1.6011 -7.755 7.27e-11 ***
sprayD -9.5833 1.6011 -5.985 9.82e-08 ***
sprayE -11.0000 1.6011 -6.870 2.75e-09 ***
sprayF 2.1667 1.6011 1.353 0.181
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.922 on 66 degrees of freedom
Multiple R-squared: 0.7244, Adjusted R-squared: 0.7036
F-statistic: 34.7 on 5 and 66 DF, p-value: < 2.2e-16
anova(m1)
2.1. Checagem de Pressupostos
# Testes Estatísticos Clássicos# H0 do Shapiro: Os dados vêm de uma distribuição Normal.shapiro.test(m1$residuals)
Shapiro-Wilk normality test
data: m1$residuals
W = 0.96006, p-value = 0.02226
# p-value 0.02 (<0.05) rejeita H0, ou seja, NÃO-NORMAL.# H0 do Bartlett: As variâncias são homogêneas (iguais) entre os grupos.bartlett.test(count ~ spray, data = insects)
Bartlett test of homogeneity of variances
data: count by spray
Bartlett's K-squared = 25.96, df = 5, p-value = 9.085e-05
# Pacote DHARMa para resíduos simulados (excelente padrão ouro moderno)library(DHARMa)
Warning: pacote 'DHARMa' foi compilado no R versão 4.5.3
This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
plot(simulateResiduals(m1)) # O gráfico aponta desvios graves, especialmente na variância.
Conclusão da Etapa 2: Os dados violam gravemente os pressupostos da ANOVA. Isso ocorre porque dados de contagem têm um limite inferior (zero) e a variância costuma crescer junto com a média. Precisamos de outras abordagens.
3. Alternativa 1: Transformação de Dados
Uma técnica clássica para lidar com dados de contagem heterocedásticos é aplicar a Transformação de Raiz Quadrada (\(\sqrt{Y}\)). Ela “comprime” os valores muito altos, estabilizando a variância.
# Criando a variável transformada (opcional se for feita direto no lm, mas bom para visualizar)insects2 <- insects |>mutate(count2 =sqrt(count)) # Ajustando novo modelo com a variável resposta dentro da raiz quadrada m2 <-lm(sqrt(count) ~ spray, data = insects)# Checando os pressupostos no novo modeloshapiro.test(m2$residuals) # Não rejeita H0 (Agora é NORMAL)
Shapiro-Wilk normality test
data: m2$residuals
W = 0.98721, p-value = 0.6814
bartlett.test(sqrt(count) ~ spray, data = insects) # Não rejeita H0 (Agora é HOMOGÊNEO)
Bartlett test of homogeneity of variances
data: sqrt(count) by spray
Bartlett's K-squared = 3.7525, df = 5, p-value = 0.5856
# Visualização da transformação (DHARMa)plot(simulateResiduals(m2)) # Sem avisos em vermelho agora! Os resíduos estão bem comportados
3.1. Comparações Múltiplas e Destransformação
Como ajustamos o modelo com sqrt(count) no R, o pacote emmeans reconhece automaticamente a transformação. Ao pedir type = "response", ele eleva as médias ao quadrado para devolver os resultados na escala original.
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
# Médias do modelo violado (Apenas para comparação pedagógica)medias <-emmeans(m1, ~spray) # Médias do modelo corrigido, na escala original destransformada medias2 <-emmeans(m2, ~spray, type ="response") # Comparações Múltiplascld(medias) # Agrupamento do modelo com pressupostos violados
cld(medias2) # Agrupamento do modelo matematicamente correto
# Nota: Destransformou os valores. O valor final aproxima-se da média aritmética, # mas não é igual devido ao efeito não-linear da raiz.# Note que houve uma ligeira diferença no agrupamento de tratamentos por significância (as letras mudaram).
4. Alternativa 2: Estatística Não-Paramétrica
Se não quisermos transformar dados, podemos usar testes não-paramétricos (baseados em postos/ranks) que não exigem normalidade, como o Kruskal-Wallis.
library(agricolae)# Teste de Kruskal-Wallisk <-kruskal(insects$count, insects$spray)k
$statistics
Chisq Df p.chisq t.value MSD
54.69134 5 1.510845e-10 1.996564 8.462804
$parameters
test p.ajusted name.t ntr alpha
Kruskal-Wallis none insects$spray 6 0.05
$means
insects.count rank std r Min Max Q25 Q50 Q75
A 14.500000 52.16667 4.719399 12 7 23 11.50 14.0 17.75
B 15.333333 54.83333 4.271115 12 7 21 12.50 16.5 17.50
C 2.083333 11.45833 1.975225 12 0 7 1.00 1.5 3.00
D 4.916667 25.58333 2.503028 12 2 12 3.75 5.0 5.00
E 3.500000 19.33333 1.732051 12 1 6 2.75 3.0 5.00
F 16.666667 55.62500 6.213378 12 9 26 12.50 15.0 22.50
$comparison
NULL
$groups
insects$count groups
F 55.62500 a
B 54.83333 a
A 52.16667 a
D 25.58333 b
E 19.33333 bc
C 11.45833 c
attr(,"class")
[1] "group"
# O agrupamento obtido aqui foi idêntico ao dos dados transformados com raiz quadrada.
A abordagem moderna e mais elegante é utilizar um GLM. Em vez de forçar os dados de contagem a se comportarem como uma curva normal (transformando-os), nós avisamos ao modelo qual é a natureza verdadeira dos dados.
Para dados de contagem, usamos a Família Poisson com a função de ligação logarítmica. \[Y \sim \text{Poisson}(\lambda)\]\[\log(\lambda) = \beta_0 + \beta_1 X_1\]
# Ajustando o modelo GLM Poisson m3 <-glm(count ~ spray, family =poisson(link ="log"), data = insects)# Apenas para ilustrar: se não especificarmos a família, o R usa a família Gaussian (Normal),# o que seria matematicamente idêntico ao nosso primeiro modelo (m1).m4 <-glm(count ~ spray, data = insects) # Checando resíduos do GLM com DHARMaplot(simulateResiduals(m3)) # Comportamento excelente!
# Obtendo médias marginais estimadas (o type = "response" aplica a exponencial para sair do log) df3 <-cld(emmeans(m3, ~ spray, type ="response"), Letters = LETTERS)
Nota: O GLM é considerado mais elegante e correto hoje em dia. Ele não adultera os dados (como na transformação) e permite o uso de inferência paramétrica poderosa. Note que o agrupamento de letras foi o mesmo da estatística não-paramétrica e do dado transformado.
5.1. Gráfico Final com Resultados do GLM
Para fechar, plotamos os resultados com o erro padrão assintótico (LCL e UCL) e as letras de significância obtidas no emmeans.
# Plotagem final corrigidadf3 |>ggplot(aes(x = spray, y = rate)) +geom_point(color ="darkred", size =3) +geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width =0.1, color ="darkred") +ylim(0, 25) +# A estética 'label' indica qual coluna contém as letras de agrupamento (.group)geom_text(aes(y = asymp.UCL, label =str_trim(.group)), vjust =-0.8, size =4) +theme_classic() +labs(x ="Inseticida",y ="Contagem de insetos vivos estimada",title ="Eficácia dos Inseticidas (GLM de Poisson)")
