Aula 6

Aula 6: ANOVA fatorial, blocos e correlação

Nesta aula trabalharemos com três situações analíticas bastante comuns em experimentação agronômica e fitopatológica:

  1. experimento fatorial, para avaliar simultaneamente os efeitos de fungicida, dose e sua interação sobre a severidade;
  2. delineamento em blocos, para analisar produtividade em experimento de campo controlando a heterogeneidade experimental;
  3. análise de associação, por meio de correlação e regressão linear simples, para investigar a relação entre variáveis quantitativas.

Ao longo do código, os comentários procuram explicar não apenas como executar os procedimentos, mas também por que eles são usados e como interpretar os resultados.

Pacotes

library(readxl)
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(emmeans)
Warning: package 'emmeans' was built under R version 4.5.3
Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'
library(multcomp)
Warning: package 'multcomp' was built under R version 4.5.3
Loading required package: mvtnorm
Warning: package 'mvtnorm' was built under R version 4.5.3
Loading required package: survival
Loading required package: TH.data
Warning: package 'TH.data' was built under R version 4.5.3
Loading required package: 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(DHARMa)
Warning: package 'DHARMa' was built under R version 4.5.3
This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library(ggthemes)
Warning: package 'ggthemes' was built under R version 4.5.3

Parte 1 — Experimento fatorial: fungicida × dose

Importação e preparo dos dados

fungicidas <- read_excel("dados_diversos.xlsx", sheet = "fungicida_vaso")

fungicidas2 <- fungicidas |> 
  mutate(
    incidencia = inf_seeds / n_seeds,
    dose = factor(dose),
    treat = factor(treat)
  )

glimpse(fungicidas2)
Rows: 20
Columns: 9
$ treat      <fct> Ionic liquid, Ionic liquid, Ionic liquid, Ionic liquid, Ion…
$ dose       <fct> 0.5, 0.5, 0.5, 0.5, 0.5, 2, 2, 2, 2, 2, 0.5, 0.5, 0.5, 0.5,…
$ rep        <dbl> 1, 2, 3, 4, 5, 2, 3, 4, 5, 1, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5
$ n_sp       <dbl> 103, 125, 210, 97, 180, 116, 166, 157, 129, 84, 121, 123, 1…
$ dis_sp     <dbl> 13, 31, 80, 28, 75, 9, 7, 12, 7, 0, 0, 6, 3, 0, 4, 3, 0, 4,…
$ n_seeds    <dbl> 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25,…
$ inf_seeds  <dbl> 10, 12, 12, 10, 11, 6, 3, 1, 7, 1, 2, 5, 1, 2, 1, 5, 2, 2, …
$ severity   <dbl> 0.12621359, 0.24800000, 0.38095238, 0.28865979, 0.41666667,…
$ incidencia <dbl> 0.40, 0.48, 0.48, 0.40, 0.44, 0.24, 0.12, 0.04, 0.28, 0.04,…

A variável incidencia foi calculada como a proporção de sementes infectadas:

\[ \text{incidência} = \frac{\text{número de sementes infectadas}}{\text{número total de sementes avaliadas}} \]

Neste primeiro conjunto, a variável resposta analisada será severity. Como treat e dose representam tratamentos experimentais, ambos devem ser tratados como fatores.

Visualização inicial

fungicidas2 |> 
  ggplot(aes(x = treat, y = severity, color = dose)) +
  geom_jitter(width = 0.08, height = 0, size = 2, alpha = 0.7) +
  facet_wrap(~ dose) +
  theme_classic() +
  labs(
    x = "Fungicida",
    y = "Severidade",
    color = "Dose"
  )

interaction.plot(
  x.factor = fungicidas2$treat,
  trace.factor = fungicidas2$dose,
  response = fungicidas2$severity,
  xlab = "Fungicida",
  ylab = "Severidade média",
  trace.label = "Dose"
)

A inspeção gráfica antecede a modelagem. Em um experimento fatorial, ela é importante porque permite avaliar se há indícios de interação entre os fatores. Em termos práticos:

  • se as linhas do gráfico de interação forem aproximadamente paralelas, a evidência visual favorece ausência de interação;
  • se as linhas se cruzarem ou divergirem claramente, isso sugere que o efeito de um fator depende do nível do outro.

Modelos candidatos

Modelo 1: efeito apenas de fungicida

m1 <- lm(severity ~ treat, data = fungicidas2)
anova(m1)
Analysis of Variance Table

Response: severity
          Df  Sum Sq Mean Sq F value   Pr(>F)   
treat      1 0.11323 0.11323  9.8892 0.005603 **
Residuals 18 0.20610 0.01145                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modelo 2: efeito apenas de dose

m2 <- lm(severity ~ dose, data = fungicidas2)
anova(m2)
Analysis of Variance Table

Response: severity
          Df   Sum Sq  Mean Sq F value  Pr(>F)  
dose       1 0.073683 0.073683  5.3991 0.03206 *
Residuals 18 0.245649 0.013647                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modelo 3: modelo fatorial com interação

m3 <- lm(severity ~ treat * dose, data = fungicidas2)
anova(m3)
Analysis of Variance Table

Response: severity
           Df   Sum Sq  Mean Sq F value    Pr(>F)    
treat       1 0.113232 0.113232  30.358 4.754e-05 ***
dose        1 0.073683 0.073683  19.755 0.0004077 ***
treat:dose  1 0.072739 0.072739  19.502 0.0004326 ***
Residuals  16 0.059678 0.003730                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

O modelo fatorial com interação pode ser escrito como:

\[ Y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \varepsilon_{ijk} \]

em que:

  • \(Y_{ijk}\) é a observação de severidade na combinação entre o fungicida \(i\) e a dose \(j\);
  • \(\mu\) é a média geral;
  • \(\alpha_i\) é o efeito do fungicida;
  • \(\beta_j\) é o efeito da dose;
  • \((\alpha\beta)_{ij}\) é o efeito de interação entre fungicida e dose;
  • \(\varepsilon_{ijk}\) é o erro aleatório, assumido com média zero e variância constante.

Na ANOVA, a estatística de teste para cada fonte de variação é dada por:

\[ F = \frac{QM_{\text{fonte}}}{QM_{\text{erro}}} \]

em que \(QM\) é o quadrado médio. Valores altos de \(F\) indicam que a variação explicada pela fonte é grande em relação à variação residual.

Interpretação da interação

Em experimentos fatoriais, a presença de interação significativa muda a forma de interpretar os efeitos principais. Quando a interação é relevante, evita-se discutir apenas “o melhor fungicida” ou “a melhor dose” isoladamente, porque o efeito de um depende do outro. Nesse caso, a interpretação correta passa a ser feita por efeitos simples, isto é, comparando doses dentro de cada fungicida, ou fungicidas dentro de cada dose.

Diagnóstico do modelo

res_m3 <- simulateResiduals(m3)
plot(res_m3)

O pacote {DHARMa} produz resíduos simulados padronizados, úteis para avaliar:

  • padrões de assimetria;
  • heterogeneidade de variâncias;
  • desvios importantes de distribuição;
  • observações discrepantes.

Embora lm() assuma erros normais com variância constante, em aplicações práticas a avaliação gráfica dos resíduos costuma ser mais informativa do que depender apenas de testes formais.

Médias ajustadas e comparação múltipla

medias_m3 <- emmeans(m3, ~ dose | treat)
medias_m3
treat = Ionic liquid:
 dose emmean     SE df lower.CL upper.CL
 0.5  0.2921 0.0273 16  0.23420   0.3500
 2    0.0501 0.0273 16 -0.00781   0.1080

treat = Tebuconazole:
 dose emmean     SE df lower.CL upper.CL
 0.5  0.0210 0.0273 16 -0.03690   0.0789
 2    0.0202 0.0273 16 -0.03768   0.0781

Confidence level used: 0.95 
cld(medias_m3, Letters = letters)
treat = Ionic liquid:
 dose emmean     SE df lower.CL upper.CL .group
 2    0.0501 0.0273 16 -0.00781   0.1080  a    
 0.5  0.2921 0.0273 16  0.23420   0.3500   b   

treat = Tebuconazole:
 dose emmean     SE df lower.CL upper.CL .group
 2    0.0202 0.0273 16 -0.03768   0.0781  a    
 0.5  0.0210 0.0273 16 -0.03690   0.0789  a    

Confidence level used: 0.95 
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.

As médias obtidas por emmeans() são médias marginais estimadas, também chamadas de médias ajustadas. Elas representam as médias previstas pelo modelo e são especialmente úteis quando o desenho não é perfeitamente balanceado ou quando há termos adicionais no modelo.

A função cld() organiza os grupos de médias por letras. Médias que compartilham a mesma letra não diferem entre si ao nível de significância adotado.

Tabela ilustrativa para a interação

Fungicida Dose 0.5 Dose 2
Ionic liquid 0.29 Aa 0.05 Ab
Tebuconazole 0.02 Ba 0.02 Aa

Leitura da tabela:

  • letras maiúsculas comparam fungicidas dentro de cada dose;
  • letras minúsculas comparam doses dentro de cada fungicida;
  • médias com a mesma letra não diferem a 5% de probabilidade.

Parte 2 — Experimento em blocos ao acaso

Importação e preparo dos dados

campo <- read_excel("dados_diversos.xlsx", sheet = "fungicida_campo")

campo2 <- campo |> 
  mutate(
    trat = factor(TRAT),
    bloco = factor(BLOCO)
  ) |> 
  transmute(
    trat,
    bloco,
    prod = PROD
  )

glimpse(campo2)
Rows: 32
Columns: 3
$ trat  <fct> 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6…
$ bloco <fct> 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2…
$ prod  <dbl> 4611.3, 4311.7, 4544.5, 3409.4, 5022.6, 4839.3, 4902.6, 4975.2, …

Neste conjunto, a variável resposta é a produtividade (prod). O fator bloco representa uma estrutura de controle experimental usada para reduzir a variabilidade não controlada entre unidades experimentais.

Visualização inicial

campo2 |> 
  ggplot(aes(x = trat, y = prod)) +
  geom_jitter(width = 0.08, alpha = 0.6, size = 2) +
  stat_summary(fun = mean, geom = "point", size = 4, shape = 18) +
  theme_classic() +
  labs(
    x = "Tratamento fungicida",
    y = expression("Produtividade (kg ha"^-1*")")
  )

Os pontos mostram os valores individuais e o losango representa a média de cada tratamento. Essa visualização ajuda a avaliar amplitude dos dados, dispersão dentro dos tratamentos e possíveis diferenças entre médias.

Modelo em blocos

m_bloco <- lm(prod ~ trat + bloco, data = campo2)
anova(m_bloco)
Analysis of Variance Table

Response: prod
          Df  Sum Sq Mean Sq F value Pr(>F)  
trat       7 2994142  427735  2.6369 0.0402 *
bloco      3  105716   35239  0.2172 0.8833  
Residuals 21 3406384  162209                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

O modelo para delineamento em blocos ao acaso pode ser escrito como:

\[ Y_{ij} = \mu + \tau_i + \beta_j + \varepsilon_{ij} \]

em que:

  • \(Y_{ij}\) é a produtividade observada no tratamento \(i\) dentro do bloco \(j\);
  • \(\mu\) é a média geral;
  • \(\tau_i\) é o efeito do tratamento;
  • \(\beta_j\) é o efeito do bloco;
  • \(\varepsilon_{ij}\) é o erro aleatório.

A principal função dos blocos é absorver parte da variabilidade sistemática do experimento. Quando os blocos capturam heterogeneidade real da área experimental, o erro residual tende a diminuir, aumentando a precisão das comparações entre tratamentos.

Diagnóstico do modelo

res_bloco <- simulateResiduals(m_bloco)
plot(res_bloco)

Médias ajustadas dos tratamentos

medias_m_bloco <- emmeans(m_bloco, ~ trat)
medias_m_bloco
 trat emmean  SE df lower.CL upper.CL
 1      4219 201 21     3800     4638
 2      4935 201 21     4516     5354
 3      5110 201 21     4691     5529
 4      5140 201 21     4722     5559
 5      5122 201 21     4703     5541
 6      5256 201 21     4838     5675
 7      5128 201 21     4709     5546
 8      5078 201 21     4659     5497

Results are averaged over the levels of: bloco 
Confidence level used: 0.95 
res_cld <- cld(medias_m_bloco, Letters = LETTERS)
res_cld
 trat emmean  SE df lower.CL upper.CL .group
 1      4219 201 21     3800     4638  A    
 2      4935 201 21     4516     5354  AB   
 8      5078 201 21     4659     5497  AB   
 3      5110 201 21     4691     5529  AB   
 5      5122 201 21     4703     5541  AB   
 7      5128 201 21     4709     5546  AB   
 4      5140 201 21     4722     5559  AB   
 6      5256 201 21     4838     5675   B   

Results are averaged over the levels of: bloco 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 8 estimates 
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.

Gráfico com médias ajustadas e grupos

res_cld |> 
  ggplot(aes(x = reorder(trat, emmean), y = emmean)) +
  geom_point(size = 3) +
  geom_text(aes(label = .group, y = emmean + 120), size = 5) +
  coord_flip() +
  theme_classic() +
  labs(
    x = "Tratamento fungicida",
    y = expression("Média ajustada de produtividade (kg ha"^-1*")")
  )

Este gráfico resume a comparação entre tratamentos já considerando o efeito de blocos no ajuste do modelo. Assim, a inferência é mais consistente do que comparar apenas médias brutas.

Parte 3 — Correlação e regressão linear

Importação dos dados

mofo <- read_excel("dados_diversos.xlsx", sheet = "mofo")

glimpse(mofo)
Rows: 52
Columns: 5
$ study <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ treat <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 1, 2, 3, 4, 5, 6, 7, …
$ inc   <dbl> 76, 53, 42, 37, 29, 42, 55, 40, 26, 18, 27, 28, 36, 76, 44, 24, …
$ scl   <dbl> 2194, 1663, 1313, 1177, 753, 1343, 1519, 516, 643, 400, 643, 921…
$ yld   <dbl> 2265, 2618, 2554, 2632, 2820, 2799, 2503, 2967, 2965, 3088, 3044…

Aqui o interesse é estudar a associação entre variáveis quantitativas. Pelo código original, a relação explorada nos gráficos é entre inc e scl. Portanto, a análise de correlação foi alinhada a essas duas variáveis para manter coerência entre teste e visualização.

Visualização geral por estudo

mofo |> 
  ggplot(aes(x = inc, y = scl, color = factor(study))) +
  geom_point(alpha = 0.7, size = 2) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_color_colorblind() +
  theme_classic() +
  labs(
    x = "Incidência",
    y = "Escala / variável scl",
    color = "Estudo"
  )
`geom_smooth()` using formula = 'y ~ x'

O coeficiente de correlação de Pearson mede a intensidade e a direção da associação linear entre duas variáveis quantitativas:

\[ r = \frac{\sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})} {\sqrt{\sum_{i=1}^{n}(x_i - \bar{x})^2} \sqrt{\sum_{i=1}^{n}(y_i - \bar{y})^2}} \]

O valor de \(r\) varia entre -1 e 1:

  • valores próximos de 1 indicam associação linear positiva forte;
  • valores próximos de -1 indicam associação linear negativa forte;
  • valores próximos de 0 indicam associação linear fraca ou ausente.

Correlação geral

cor.test(mofo$inc, mofo$scl)

    Pearson's product-moment correlation

data:  mofo$inc and mofo$scl
t = 5.4615, df = 50, p-value = 1.484e-06
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4061148 0.7577627
sample estimates:
      cor 
0.6112689

O teste associado avalia a hipótese nula:

\[ H_0: \rho = 0 \]

contra a hipótese alternativa de que a correlação populacional difere de zero.

Correlações por estudo

correlacoes <- mofo |> 
  group_by(study) |> 
  summarise(r = cor(inc, scl, use = "complete.obs"))

correlacoes
# A tibble: 4 × 2
  study     r
  <dbl> <dbl>
1     1 0.913
2     2 0.849
3     3 0.844
4     4 0.927

Inspeção das distribuições

hist(mofo$inc, main = "Histograma de incidência", xlab = "inc")

hist(mofo$scl, main = "Histograma de scl", xlab = "scl")

A inspeção visual das distribuições é útil porque correlação de Pearson e regressão linear simples assumem, em termos de interpretação inferencial, relação aproximadamente linear e ausência de distorções severas causadas por valores extremos.

Gráfico final com coeficientes por estudo

label_data <- mofo |>
  group_by(study) |>
  summarise(
    x_pos = max(inc, na.rm = TRUE),
    y_pos = predict(
      lm(scl ~ inc, data = cur_data()),
      newdata = data.frame(inc = max(inc, na.rm = TRUE))
    ) |> as.numeric(),
    .groups = "drop"
  ) |>
  left_join(correlacoes, by = "study")
Warning: There was 1 warning in `summarise()`.
ℹ In argument: `y_pos = as.numeric(...)`.
ℹ In group 1: `study = 1`.
Caused by warning:
! `cur_data()` was deprecated in dplyr 1.1.0.
ℹ Please use `pick()` instead.
mofo |> 
  ggplot(aes(x = inc, y = scl, color = factor(study))) +
  geom_point(alpha = 0.5, size = 2) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 1) +
  geom_text(
    data = label_data,
    aes(x = x_pos, y = y_pos, label = paste0("r = ", format(round(r, 2), nsmall = 2))),
    hjust = -0.15,
    vjust = 0.5,
    fontface = "italic",
    show.legend = FALSE
  ) +
  scale_color_colorblind() +
  scale_x_continuous(expand = expansion(mult = c(0.05, 0.20))) +
  theme_classic(base_size = 12) +
  theme(
    legend.position = "top",
    axis.title = element_text(face = "bold"),
    axis.text = element_text(color = "black")
  ) +
  labs(
    x = "Incidência da doença (%)",
    y = "scl",
    color = "Estudo"
  )
`geom_smooth()` using formula = 'y ~ x'

Esse gráfico reúne três elementos importantes:

  1. dados observados, representados pelos pontos;
  2. tendência linear ajustada, representada pelas retas;
  3. coeficiente de correlação por estudo, indicado ao final de cada reta.

Essa abordagem é útil quando se deseja mostrar que a associação pode variar entre estudos, experimentos ou ambientes.

Síntese conceitual da aula

Nesta aula, três ideias estatísticas centrais foram trabalhadas.

1. Interação em experimentos fatoriais

Quando há interação, o efeito de um fator depende do nível do outro. Nessa situação, a interpretação deve se concentrar em comparações dentro de cada combinação relevante.

2. Uso de blocos para ganhar precisão

O bloco não é o foco biológico principal, mas uma ferramenta de controle experimental. Seu papel é reduzir a variabilidade residual e melhorar a sensibilidade do teste para detectar diferenças entre tratamentos.

3. Correlação não implica causalidade

Mesmo quando a correlação é alta, isso indica associação estatística, e não necessariamente uma relação causal. A interpretação deve sempre considerar o contexto biológico e experimental.