Análise Estatística: Dados de Contagem de Insetos

# Introdução

#Este documento apresenta o fluxo de análise estatística para dados de contagem, utilizando o dataset `InsectSprays`. O objetivo é comparar a eficácia de diferentes sprays inseticidas, passando pela verificação de pressupostos e modelos lineares generalizados.

# 1. Preparação e Exploração

#Primeiro, carregamos as bibliotecas necessárias para diagnóstico, modelagem e visualização.

#| message: false
#| warning: false

library(tidyverse)    # Manipulação e ggplot2

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.6.0
✔ ggplot2   4.0.1     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.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

library(performance)  # Diagnóstico de modelos (AIC, R2, pressupostos)
library(DHARMa)       # Resíduos simulados para modelos não-normais

This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')

library(emmeans)      # Médias estimadas e comparações

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

library(multcomp)     # Testes de comparação múltipla (CLD)

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

Attaching package: 'MASS'

The following object is masked from 'package:dplyr':

    select


Attaching package: 'TH.data'

The following object is masked from 'package:MASS':

    geyser

library(multcompView) # Suporte para letras de compressão
library(agricolae)    # Testes não-paramétricos (Kruskal-Wallis)

# Carregando dados nativos do R
insetos <- InsectSprays

0.0.1 Visualização Inicial

A análise de contagem geralmente revela que grupos com médias maiores também possuem maior variabilidade (heterocedasticidade).

ggplot(insetos, aes(x = spray, y = count, fill = spray)) +
  geom_boxplot(outlier.color = NA, alpha = 0.5) +
  geom_jitter(width = 0.1, alpha = 0.6) +
  theme_minimal() +
  labs(title = "Contagem de Insetos por Tipo de Spray",
       y = "Número de Insetos Vivos", x = "Tipo de Spray")

1 2. Modelo Linear Clássico (ANOVA)

Tentamos ajustar um modelo linear simples (\(y = \mu + \tau_i + \epsilon\)).

m1 <- lm(count ~ spray, data = insetos)

# Checagem de Pressupostos
# Em dados de contagem, a variância costuma aumentar com a média (Poisson)
check_normality(m1)

Warning: Non-normality of residuals detected (p = 0.022).

check_heteroscedasticity(m1)

Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).

Se o teste de Breusch-Pagan ou Levene indicar heterocedasticidade, o valor de p da ANOVA pode ser espúrio.

2 3. Tratamento de Dados: Transformação

A transformação raiz quadrada (\(\sqrt{x}\)) é a técnica clássica para “estabilizar” a variância em dados de Poisson.

m2 <- lm(sqrt(count) ~ spray, data = insetos)

# Verificando se a transformação corrigiu o problema
check_heteroscedasticity(m2) 

OK: Error variance appears to be homoscedastic (p = 0.854).

check_normality(m2)

OK: residuals appear as normally distributed (p = 0.681).

# ANOVA do modelo transformado
anova(m2)

Analysis of Variance Table

Response: sqrt(count)
          Df Sum Sq Mean Sq F value    Pr(>F)    
spray      5 88.438 17.6876  44.799 < 2.2e-16 ***
Residuals 66 26.058  0.3948                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3 4. Abordagem Generalizada (GLM Poisson)

A alternativa moderna é não transformar os dados, mas usar um modelo que entenda a distribuição de Poisson e utilize uma função de ligação log.

# Ajuste do GLM (Família Poisson, Link Log)
m3 <- glm(count ~ spray, 
          family = poisson(link = "log"), 
          data = insetos)

# Diagnóstico com resíduos simulados (DHARMa)
# Essencial para detectar Overdispersion (variância > média)
res_m3 <- simulateResiduals(m3)
plot(res_m3)

4 5. Médias Estimadas e Comparação Múltipla

Para reportar os resultados do GLM, destransformamos os valores da escala logarítmica de volta para a escala original de contagem usando type = "response".

# Obtendo médias destransformadas e letras do teste de Tukey
medias_glm <- emmeans(m3, ~ spray, type = "response")
cld_glm <- cld(medias_glm, Letters = LETTERS, adjust = "tukey")

# Exibição da tabela de resultados
print(cld_glm)

 spray  rate    SE  df asymp.LCL asymp.UCL .group
 C      2.08 0.417 Inf      1.23      3.53  A    
 E      3.50 0.540 Inf      2.33      5.25  AB   
 D      4.92 0.640 Inf      3.49      6.93   B   
 A     14.50 1.100 Inf     11.88     17.70    C  
 B     15.33 1.130 Inf     12.63     18.62    C  
 F     16.67 1.180 Inf     13.84     20.07    C  

Confidence level used: 0.95 
Conf-level adjustment: sidak method for 6 estimates 
Intervals are back-transformed from the log scale 
P value adjustment: tukey method for comparing a family of 6 estimates 
Tests are performed on the log scale 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

5 6. Visualização dos Resultados Finais

Gráfico de pontos com intervalos de confiança (95%) e as letras que indicam diferenças significativas (letras diferentes = diferença estatística).

ggplot(cld_glm, aes(x = reorder(spray, rate), y = rate)) +
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.1, color = "gray30") +
  geom_point(size = 4, color = "darkblue") +
  geom_text(aes(label = .group, y = asymp.UCL), vjust = -1.2, size = 4, fontface = "bold") +
  scale_y_continuous(limits = c(0, 25)) +
  theme_classic() +
  labs(x = "Tipo de Spray",
       y = "Média de Insetos (Escala Real)",
       title = "Eficácia Comparativa dos Inseticidas",
       subtitle = "Estimativas obtidas via GLM Poisson (Link Log)") +
  coord_flip()