ANOVA e Planejamento de Experimentos (PBL)

.title[
# ANOVA e Planejamento de Experimentos (PBL)
]
.subtitle[
## Aplicações em Saúde e Agronomia — teoria + prática em R
]
.author[
### Prof. Hidelbrando Ferreira Rodrigues
]
.date[
### 27/04/2026
]

---

# ANOVA e Planejamento de Experimentos (PBL)
## Saúde e Agronomia: mesmo modelo, diferentes unidades experimentais

---

## Objetivos da aula

- Entender a lógica da **ANOVA** como comparação **global** de médias.
- Conectar ANOVA ao **planejamento**: randomização, replicação e blocagem.
- Aplicar **delineamentos clássicos**:
  - DIC, DBC, Quadrado Latino.
- Tomar decisões pós-ANOVA:
  - fatores **qualitativos** → testes de médias / contrastes,
  - fatores **quantitativos** → **regressão** (com teste de desvios).

---

## Estrutura PBL

| Etapa | O que fazemos |
|------:|---------------|
| **1. Problema** | caso realista — contexto e pergunta |
| **2. Delineamento** | hipóteses, UE, fator, blocos |
| **3. Modelo + ANOVA** | F global e tabela de ANOVA |
| **4. Diagnóstico** | verificação de pressupostos |
| **5. Pós-ANOVA** | médias / contrastes / regressão |
| **6. Recomendação** | conclusão aplicada |

---

# Parte 1 — ANOVA: ideia e mecânica

---

## Motivação — Saúde: 3 protocolos, 4 hospitais

- Protocolos A, B, C aplicados em pacientes
- Desfecho contínuo: escore funcional após 30 dias
- Hospitais com perfis distintos (infraestrutura, gravidade dos casos)

---

## Motivação — Agronomia: 4 épocas de corte + gradiente

- Tratamentos A, B, C, D (épocas ou cultivares)
- Gradiente de solo/declividade no campo experimental
- Blocos para absorver esse gradiente e reduzir `$QM_E$`

---

## ANOVA one-way — Modelo

`$$Y_{ij} = \mu + \tau_i + \varepsilon_{ij}$$`

- `$i = 1,\dots,k$`: tratamentos; `$\;\; j = 1,\dots,n_i$`: repetições
- `$\varepsilon_{ij} \sim N(0,\,\sigma^2)$` independentes
- `$H_0: \mu_1 = \mu_2 = \cdots = \mu_k$`
- `$F = \dfrac{QM_{Trat}}{QM_E} \sim F_{k-1,\; N-k}$` sob `$H_0$`

---

## Como ler a tabela ANOVA

| Coluna | O que representa | Cálculo |
|--------|-----------------|---------|
| **FV** | Fonte de Variação | Tratamentos / Resíduo / Total |
| **GL** | Graus de liberdade | `$k-1$` / `$N-k$` / `$N-1$` |
| **SQ** | Soma de Quadrados | variação atribuída a cada FV |
| **QM** | Quadrado Médio | `$SQ \div GL$` |
| **F**  | Estatística de teste | `$QM_{Trat} \div QM_{Res}$` |
| **p**  | Valor-p | `$P(F \geq F_{obs} \mid H_0)$` |

.box-azul[
💡 **Regra prática:** p < 0,05 → rejeita-se `$H_0$` → ao menos um par de médias difere.
Um F grande indica que a variação *entre* tratamentos supera a variação *dentro* deles.
]

---

## Pressupostos da ANOVA

1. **Independência** — garantida pela randomização
2. **Normalidade** dos resíduos
3. **Homogeneidade** de variâncias

**Checagem:**
- QQ-plot e resíduos vs ajustados (gráficos — sempre priorize)
- Shapiro-Wilk · Levene / Bartlett (testes formais)

.box-laranja[
⚠️ **Atenção — Shapiro-Wilk com n grande:** com amostras grandes o teste tem poder
muito elevado e detecta desvios triviais. Priorize sempre a análise gráfica.
]

---

# Parte 2 — PBL 1 (Saúde): DIC vs DBC

---

## Caso PBL 1 — Protocolo Hospitalar

3 protocolos (A, B, C) aplicados em pacientes de 4 hospitais.  
Desfecho: **escore funcional** aos 30 dias.

**Tarefas:**
- Definir unidade experimental, fator e níveis.
- Decidir: hospital entra como **bloco**?
- Ajustar DIC e DBC no R e comparar `$QM_{Res}$`.

---

## Dados simulados (Saúde)

``` r
set.seed(10)

dados_saude <- expand.grid(
  protocolo = c("A", "B", "C"),
  hospital  = paste0("H", 1:4),
  rep       = 1:12
)

efeito_prot <- c(A = 0, B = 4, C = 7)
efeito_hosp <- c(H1 = -3, H2 = 0, H3 = 2, H4 = 5)

dados_saude$escore <- 60 +
  efeito_prot[dados_saude$protocolo] +
  efeito_hosp[dados_saude$hospital] +
  rnorm(nrow(dados_saude), 0, 6)

dados_saude$protocolo <- factor(dados_saude$protocolo)
dados_saude$hospital  <- factor(dados_saude$hospital)
```
]

``` r
dplyr::glimpse(dados_saude)
```

```
## Rows: 144
## Columns: 4
## $ protocolo <fct> A, B, C, A, B, C, A, B, C, A, B, C, A, B, C, A, B, C, A, B, …
## $ hospital  <fct> H1, H1, H1, H2, H2, H2, H3, H3, H3, H4, H4, H4, H1, H1, H1, …
## $ rep       <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, …
## $ escore    <dbl> 57.11248, 59.89448, 55.77202, 56.40499, 65.76727, 69.33877, …
```

]

---

## DIC vs DBC (Saúde)

``` r
m_dic_saude <- aov(escore ~ protocolo,            data = dados_saude)
m_dbc_saude <- aov(escore ~ protocolo + hospital, data = dados_saude)

summary(m_dic_saude)
```

```
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## protocolo     2   1384   691.9   16.12 4.96e-07 ***
## Residuals   141   6052    42.9                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

``` r
summary(m_dbc_saude)
```

```
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## protocolo     2   1384   691.9   21.38 8.14e-09 ***
## hospital      3   1587   529.0   16.35 3.80e-09 ***
## Residuals   138   4465    32.4                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

.box-azul[
💡 Observe como o `$QM_{Res}$` diminui no DBC ao absorver a variação entre hospitais,
aumentando o poder do teste F sobre os protocolos.
]

---

## Diagnóstico (Saúde)

``` r
plot_diagnosticos(m_dbc_saude, "Diagnóstico — DBC Saúde")
```

![](Análise-de-Variância-e-Planejamento-de-Experimentos_files/figure-html/saude-diagnostico-1.png)

``` r
shapiro.test(residuals(m_dbc_saude))   # n = 144: priorize o QQ-plot
```

```
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(m_dbc_saude)
## W = 0.98104, p-value = 0.04355
```

``` r
car::leveneTest(escore ~ protocolo, data = dados_saude)
```

```
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   2  0.4038 0.6685
##       141
```

---

# Parte 3 — DIC: ANOVA "na mão" e no R (peixes)

---

## Dados (peixes)

``` r
dados_peixe <- tibble::tribble(
  ~racao, ~peso,
  "A", 1.04,  "A", 1.33,  "A", 1.14,
  "B", 1.12,  "B", 0.98,  "B", 1.14,
  "C", 1.40,  "C", 1.22,  "C", 1.24,
  "D", 0.93,  "D", 1.21,  "D", 1.16
) |>
  mutate(racao = factor(racao))

dados_peixe
```

```
## # A tibble: 12 × 2
##    racao  peso
##    <fct> <dbl>
##  1 A      1.04
##  2 A      1.33
##  3 A      1.14
##  4 B      1.12
##  5 B      0.98
##  6 B      1.14
##  7 C      1.4 
##  8 C      1.22
##  9 C      1.24
## 10 D      0.93
## 11 D      1.21
## 12 D      1.16
```

---

## ANOVA DIC "na mão"

``` r
anova_dic_manual <- function(y, trat) {
  trat <- factor(trat); N <- length(y); k <- nlevels(trat)
  n_i  <- as.numeric(table(trat))
  if (length(unique(n_i)) != 1) stop("Função assume DIC balanceado.")
  tot <- tapply(y, trat, sum); G <- sum(y); C <- G^2 / N; n <- n_i[1]
  SQtot  <- sum(y^2) - C
  SQtrat <- sum(tot^2) / n - C
  SQres  <- SQtot - SQtrat
  gl_trat <- k - 1; gl_res <- N - k
  QMtrat  <- SQtrat / gl_trat; QMres <- SQres / gl_res
  Fcalc   <- QMtrat / QMres
  pval    <- pf(Fcalc, gl_trat, gl_res, lower.tail = FALSE)
  CV      <- sqrt(QMres) / mean(y) * 100
  list(
    tabela = tibble::tibble(
      FV = c("Tratamentos", "Resíduo", "Total"),
      GL = c(gl_trat, gl_res, N - 1),
      SQ = round(c(SQtrat, SQres, SQtot), 4),
      QM = round(c(QMtrat, QMres, NA_real_), 4),
      F  = round(c(Fcalc,  NA_real_, NA_real_), 2),
      p  = round(c(pval,   NA_real_, NA_real_), 4)
    ), CV = round(CV, 2)
  )
}
res <- anova_dic_manual(dados_peixe$peso, dados_peixe$racao)
res$tabela; cat("CV (%) =", res$CV)
```

```
## # A tibble: 3 × 6
##   FV             GL     SQ      QM     F      p
##   <chr>       <dbl>  <dbl>   <dbl> <dbl>  <dbl>
## 1 Tratamentos     3 0.0784  0.0261   1.7  0.243
## 2 Resíduo         8 0.123   0.0153  NA   NA    
## 3 Total          11 0.201  NA       NA   NA
```

```
## CV (%) = 10.68
```

---

## ANOVA com `aov()` + CV

``` r
m_peixe <- aov(peso ~ racao, data = dados_peixe)
summary(m_peixe)
```

```
##             Df  Sum Sq Mean Sq F value Pr(>F)
## racao        3 0.07842 0.02614   1.705  0.243
## Residuals    8 0.12267 0.01533
```

``` r
tab   <- summary(m_peixe)[[1]]
qmres <- tab["Residuals", "Mean Sq"]
cv    <- sqrt(qmres) / mean(dados_peixe$peso) * 100
cat("CV (%) =", round(cv, 2))
```

```
## CV (%) = 10.68
```

---

# Parte 4 — DBC: épocas de corte de alfafa

---

## Dados (alfafa)

``` r
dados_alfafa <- tibble::tribble(
  ~TRAT, ~BLOCO, ~PROD,
  "B","I",1.58, "D","I",2.56, "C","I",2.29, "A","I",2.89,
  "C","II",2.98,"A","II",2.88,"D","II",2.00,"B","II",1.28,
  "B","III",1.22,"C","III",1.55,"A","III",1.88,"D","III",1.82,
  "A","IV",2.90,"D","IV",2.20,"C","IV",1.95,"B","IV",1.21,
  "C","V",1.15, "B","V",1.30, "D","V",1.33, "A","V",2.20,
  "D","VI",1.00,"A","VI",2.65,"B","VI",1.66,"C","VI",1.12
) |>
  mutate(
    TRAT  = factor(TRAT),
    BLOCO = factor(BLOCO, levels = c("I","II","III","IV","V","VI"))
  )

dados_alfafa
```

```
## # A tibble: 24 × 3
##    TRAT  BLOCO  PROD
##    <fct> <fct> <dbl>
##  1 B     I      1.58
##  2 D     I      2.56
##  3 C     I      2.29
##  4 A     I      2.89
##  5 C     II     2.98
##  6 A     II     2.88
##  7 D     II     2   
##  8 B     II     1.28
##  9 B     III    1.22
## 10 C     III    1.55
## # ℹ 14 more rows
```

---

## ANOVA DBC + CV + Tukey (alfafa)

``` r
m_alfafa <- aov(PROD ~ TRAT + BLOCO, data = dados_alfafa)
summary(m_alfafa)
```

```
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## TRAT         3  4.382  1.4607   8.471 0.00157 **
## BLOCO        5  2.759  0.5518   3.200 0.03656 * 
## Residuals   15  2.587  0.1724                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

``` r
tab_alf   <- summary(m_alfafa)[[1]]
qmres_alf <- tab_alf["Residuals", "Mean Sq"]
glres_alf <- tab_alf["Residuals", "Df"]
cv        <- sqrt(qmres_alf) / mean(dados_alfafa$PROD) * 100
cat("CV (%) =", round(cv, 2))
```

```
## CV (%) = 21.86
```

``` r
tuk_alf <- agricolae::HSD.test(m_alfafa, "TRAT", alpha = 0.05, group = TRUE)
tuk_alf$groups
```

```
##       PROD groups
## A 2.566667      a
## C 1.840000      b
## D 1.818333      b
## B 1.375000      b
```

---

## Diagnóstico (alfafa)

``` r
plot_diagnosticos(m_alfafa, "Diagnóstico — DBC Alfafa")
```

![](Análise-de-Variância-e-Planejamento-de-Experimentos_files/figure-html/dbc-alfafa-diagnostico-1.png)

``` r
shapiro.test(residuals(m_alfafa))
```

```
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(m_alfafa)
## W = 0.97524, p-value = 0.7948
```

``` r
car::leveneTest(PROD ~ TRAT, data = dados_alfafa)
```

```
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3   2.056 0.1384
##       20
```

---

## Médias com IC 95% (alfafa)

``` r
r_alf <- nlevels(dados_alfafa$BLOCO)

dados_alfafa |>
  group_by(TRAT) |>
  summarise(media = mean(PROD), .groups = "drop") |>
  mutate(
    se   = sqrt(qmres_alf / r_alf),          # EP da média no DBC
    ic95 = qt(0.975, glres_alf) * se
  ) |>
  ggplot(aes(TRAT, media)) +
  geom_col(fill = "#1B6CA8") +
  geom_errorbar(aes(ymin = media - ic95, ymax = media + ic95), width = 0.25) +
  geom_text(aes(label = round(media, 2)), vjust = -1.3) +
  labs(
    title    = "DBC — Alfafa: produtividade média por época de corte",
    subtitle = "Barras de erro: IC 95% baseado no QM do resíduo",
    x = "Tratamento", y = "Produtividade (kg/parcela)"
  ) +
  theme_minimal(base_size = 16)
```

![](Análise-de-Variância-e-Planejamento-de-Experimentos_files/figure-html/dbc-alfafa-plot-1.png)

---

# Parte 5 — DBC: diâmetro de mudas de laranjeira

---

## Dados (laranjeira)

``` r
dados_laranja <- tibble::tribble(
  ~TRAT, ~BLOCO, ~VR,
  "T1","I",1.75, "T1","II",2.03, "T1","III",2.12, "T1","IV",2.14,
  "T2","I",2.05, "T2","II",2.26, "T2","III",2.42, "T2","IV",2.53,
  "T3","I",2.34, "T3","II",2.02, "T3","III",2.43, "T3","IV",2.26,
  "T4","I",2.80, "T4","II",3.84, "T4","III",3.44, "T4","IV",3.09,
  "T5","I",1.95, "T5","II",2.15, "T5","III",1.99, "T5","IV",2.17,
  "T6","I",3.51, "T6","II",3.32, "T6","III",3.68, "T6","IV",3.31
) |>
  mutate(
    TRAT  = factor(TRAT),
    BLOCO = factor(BLOCO, levels = c("I","II","III","IV"))
  )

dados_laranja
```

```
## # A tibble: 24 × 3
##    TRAT  BLOCO    VR
##    <fct> <fct> <dbl>
##  1 T1    I      1.75
##  2 T1    II     2.03
##  3 T1    III    2.12
##  4 T1    IV     2.14
##  5 T2    I      2.05
##  6 T2    II     2.26
##  7 T2    III    2.42
##  8 T2    IV     2.53
##  9 T3    I      2.34
## 10 T3    II     2.02
## # ℹ 14 more rows
```

---

## ANOVA DBC + Tukey (laranjeira)

``` r
m_laranja <- aov(VR ~ TRAT + BLOCO, data = dados_laranja)
summary(m_laranja)
```

```
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## TRAT         5  8.133  1.6267   30.43 2.42e-07 ***
## BLOCO        3  0.253  0.0845    1.58    0.236    
## Residuals   15  0.802  0.0535                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

``` r
tab_lar   <- summary(m_laranja)[[1]]
qmres_lar <- tab_lar["Residuals", "Mean Sq"]
glres_lar <- tab_lar["Residuals", "Df"]

tuk_lar <- agricolae::HSD.test(m_laranja, "TRAT", alpha = 0.05, group = TRUE)
tuk_lar$groups
```

```
##        VR groups
## T6 3.4550      a
## T4 3.2925      a
## T2 2.3150      b
## T3 2.2625      b
## T5 2.0650      b
## T1 2.0100      b
```

---

## Diagnóstico (laranjeira)

``` r
plot_diagnosticos(m_laranja, "Diagnóstico — DBC Laranjeira")
```

![](Análise-de-Variância-e-Planejamento-de-Experimentos_files/figure-html/dbc-laranja-diagnostico-1.png)

``` r
shapiro.test(residuals(m_laranja))
```

```
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(m_laranja)
## W = 0.96675, p-value = 0.5879
```

``` r
car::leveneTest(VR ~ TRAT, data = dados_laranja)
```

```
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  5   2.256 0.09289 .
##       18                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---

## Médias com IC 95% (laranjeira)

``` r
r_lar <- nlevels(dados_laranja$BLOCO)

dados_laranja |>
  group_by(TRAT) |>
  summarise(media = mean(VR), .groups = "drop") |>
  mutate(
    se   = sqrt(qmres_lar / r_lar),
    ic95 = qt(0.975, glres_lar) * se
  ) |>
  ggplot(aes(TRAT, media)) +
  geom_col(fill = "#1B6CA8") +
  geom_errorbar(aes(ymin = media - ic95, ymax = media + ic95), width = 0.25) +
  geom_text(aes(label = round(media, 2)), vjust = -1.3) +
  labs(
    title    = "DBC — Laranjeira: diâmetro médio por tratamento",
    subtitle = "Barras de erro: IC 95% baseado no QM do resíduo",
    x = "Tratamento", y = "Diâmetro (mm)"
  ) +
  theme_minimal(base_size = 16)
```

![](Análise-de-Variância-e-Planejamento-de-Experimentos_files/figure-html/dbc-laranja-plot-1.png)

---

# Parte 6 — Quadrado Latino (suínos)

---

## Dados (suíno)

``` r
dados_suino <- tibble::tribble(
  ~LIN,  ~COL,  ~TRAT, ~VR,
  "R1","F1","A",35, "R1","F2","B",33, "R1","F3","D",28, "R1","F4","C",28,
  "R2","F1","B",15, "R2","F2","C",40, "R2","F3","A",29, "R2","F4","D",14,
  "R3","F1","C",31, "R3","F2","D",36, "R3","F3","B",20, "R3","F4","A",27,
  "R4","F1","D",19, "R4","F2","A",46, "R4","F3","C",39, "R4","F4","B",12
) |>
  mutate(across(c(LIN, COL, TRAT), factor))

dados_suino
```

```
## # A tibble: 16 × 4
##    LIN   COL   TRAT     VR
##    <fct> <fct> <fct> <dbl>
##  1 R1    F1    A        35
##  2 R1    F2    B        33
##  3 R1    F3    D        28
##  4 R1    F4    C        28
##  5 R2    F1    B        15
##  6 R2    F2    C        40
##  7 R2    F3    A        29
##  8 R2    F4    D        14
##  9 R3    F1    C        31
## 10 R3    F2    D        36
## 11 R3    F3    B        20
## 12 R3    F4    A        27
## 13 R4    F1    D        19
## 14 R4    F2    A        46
## 15 R4    F3    C        39
## 16 R4    F4    B        12
```

---

## ANOVA (Quadrado Latino) + Tukey

``` r
m_suino <- aov(VR ~ TRAT + LIN + COL, data = dados_suino)
summary(m_suino)
```

```
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## TRAT         3  636.5  212.17  35.361 0.000329 ***
## LIN          3   89.0   29.67   4.944 0.046240 *  
## COL          3  741.5  247.17  41.194 0.000213 ***
## Residuals    6   36.0    6.00                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

``` r
tab_sui   <- summary(m_suino)[[1]]
qmres_sui <- tab_sui["Residuals", "Mean Sq"]
glres_sui <- tab_sui["Residuals", "Df"]
cv        <- sqrt(qmres_sui) / mean(dados_suino$VR) * 100
cat("CV (%) =", round(cv, 2))
```

```
## CV (%) = 8.67
```

``` r
tuk_sui <- agricolae::HSD.test(m_suino, "TRAT", alpha = 0.05, group = TRUE)
tuk_sui$groups
```

```
##      VR groups
## C 34.50      a
## A 34.25      a
## D 24.25      b
## B 20.00      b
```

---

## Diagnóstico (suíno)

``` r
plot_diagnosticos(m_suino, "Diagnóstico — Quadrado Latino (Suíno)")
```

![](Análise-de-Variância-e-Planejamento-de-Experimentos_files/figure-html/dql-suino-diagnostico-1.png)

``` r
shapiro.test(residuals(m_suino))
```

```
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(m_suino)
## W = 0.98944, p-value = 0.9989
```

``` r
car::leveneTest(VR ~ TRAT, data = dados_suino)
```

```
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.2445 0.8636
##       12
```

---

## Médias com IC 95% (suíno)

``` r
p_sui <- nlevels(dados_suino$LIN)  # No DQL p×p: cada trat. aparece p vezes

dados_suino |>
  group_by(TRAT) |>
  summarise(media = mean(VR), .groups = "drop") |>
  mutate(
    se   = sqrt(qmres_sui / p_sui),
    ic95 = qt(0.975, glres_sui) * se
  ) |>
  ggplot(aes(TRAT, media)) +
  geom_col(fill = "#1B6CA8") +
  geom_errorbar(aes(ymin = media - ic95, ymax = media + ic95), width = 0.25) +
  geom_text(aes(label = round(media, 1)), vjust = -1.3) +
  labs(
    title    = "DQL — Suíno: ganho médio por tratamento (ração)",
    subtitle = "Barras de erro: IC 95% baseado no QM do resíduo",
    x = "Tratamento (ração)", y = "Ganho de peso (g)"
  ) +
  theme_minimal(base_size = 16)
```

![](Análise-de-Variância-e-Planejamento-de-Experimentos_files/figure-html/dql-suino-plot-1.png)

---

# Parte 7 — Pós-ANOVA: testes de médias, contrastes e Dunnett (IVA)

---

## Dados (IVA)

``` r
dados_iva <- tibble::tribble(
  ~TRAT, ~VR,
  "A",40.2, "A",49.3, "A",40.1, "A",43.0, "A",52.4,
  "B",42.0, "B",44.5, "B",53.0, "B",54.5, "B",51.0,
  "C",47.1, "C",55.5, "C",58.3, "C",53.4, "C",45.7,
  "D",38.1, "D",45.9, "D",43.7, "D",40.6, "D",36.7
) |>
  mutate(TRAT = factor(TRAT))
```

---

## ANOVA + diagnóstico (IVA)

``` r
m_iva <- aov(VR ~ TRAT, data = dados_iva)
anova(m_iva)
```

```
## Analysis of Variance Table
## 
## Response: VR
##           Df Sum Sq Mean Sq F value Pr(>F)  
## TRAT       3 343.75 114.583  4.3676 0.0199 *
## Residuals 16 419.76  26.235                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

``` r
plot_diagnosticos(m_iva, "Diagnóstico — DIC IVA")
```

![](Análise-de-Variância-e-Planejamento-de-Experimentos_files/figure-html/dic-iva-diagnostico-1.png)

``` r
shapiro.test(residuals(m_iva))
```

```
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(m_iva)
## W = 0.91762, p-value = 0.08919
```

``` r
car::leveneTest(VR ~ TRAT, data = dados_iva)
```

```
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.2031 0.8927
##       16
```

---

## Tukey / SNK / LSD (IVA)

``` r
tab   <- anova(m_iva)
glres <- tab$Df[2];  qmres <- tab$`Mean Sq`[2]

tuk <- agricolae::HSD.test(m_iva, "TRAT", alpha = 0.05, group = TRUE)
tuk$groups
```

```
##   VR groups
## C 52      a
## B 49     ab
## A 45     ab
## D 41      b
```

``` r
snk <- agricolae::SNK.test(
  y = dados_iva$VR, trt = dados_iva$TRAT,
  DFerror = glres, MSerror = qmres, alpha = 0.05, group = TRUE
)
snk$groups
```

```
##   dados_iva$VR groups
## C           52      a
## B           49     ab
## A           45     ab
## D           41      b
```

``` r
lsd <- agricolae::LSD.test(m_iva, "TRAT", alpha = 0.05, group = TRUE)
lsd$groups
```

```
##   VR groups
## C 52      a
## B 49     ab
## A 45     bc
## D 41      c
```

.box-laranja[
⚠️ **LSD sem proteção:** não controla a taxa global de erro tipo I (FWER).
Use apenas quando F global for significativo **e** o número de comparações for pequeno (≤ 3 pares).
]

---

## Scott-Knott (IVA)

``` r
if (requireNamespace("ScottKnott", quietly = TRUE)) {
  sk_iva <- ScottKnott::SK(m_iva, which = "TRAT")
  summary(sk_iva)
} else {
  message("Instale o pacote: install.packages('ScottKnott')")
}
```

```
## List of 2
##  $ new.fact.lev:List of 1
##   ..$ TRAT: chr [1:4] "A" "B" "C" "D"
##  $ grid.data   :'data.frame':    4 obs. of  1 variable:
##   ..$ TRAT: chr [1:4] "A" "B" "C" "D"
## Goups of means at sig.level = 0.05 
##   Means G1 G2
## C 52.00  a   
## B 49.00  a   
## A 45.00     b
## D 41.00     b
```

---

## Contrastes planejados (IVA)

``` r
dados_iva2 <- dados_iva

Cmat <- matrix(c(
   1/3, -1,   1/3,  1/3,   # C1: (A+C+D)/3 vs B
   1/2,  0,  -1,    1/2,   # C2: (A+D)/2 vs C
   1,    0,   0,   -1      # C3: A vs D
), nrow = 4, ncol = 3, byrow = FALSE)

rownames(Cmat) <- levels(dados_iva2$TRAT)
colnames(Cmat) <- c("C1", "C2", "C3")

# Ortogonalidade: fora da diagonal deve ser zero
round(crossprod(Cmat), 6)
```

```
##          C1  C2 C3
## C1 1.333333 0.0  0
## C2 0.000000 1.5  0
## C3 0.000000 0.0  2
```

]

``` r
contrasts(dados_iva2$TRAT) <- Cmat
m_iva_con <- aov(VR ~ TRAT, data = dados_iva2)
summary(m_iva_con,
        split = list(TRAT = list(C1 = 1, C2 = 2, C3 = 3)))
```

```
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## TRAT         3  343.8  114.58   4.368 0.01990 * 
##   TRAT: C1   1   33.8   33.75   1.286 0.27341   
##   TRAT: C2   1  270.0  270.00  10.292 0.00548 **
##   TRAT: C3   1   40.0   40.00   1.525 0.23474   
## Residuals   16  419.8   26.24                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

.box-verde[
✅ Contrastes ortogonais particionam `$SQ_{Trat}$` em componentes independentes —
cada contraste responde a uma pergunta distinta.
]

]

---

## Dunnett (IVA) — vs controle A

``` r
dados_iva3      <- dados_iva
dados_iva3$TRAT <- relevel(dados_iva3$TRAT, ref = "A")

m_iva_dun <- aov(VR ~ TRAT, data = dados_iva3)

dun <- multcomp::glht(
  m_iva_dun,
  linfct = multcomp::mcp(TRAT = "Dunnett")
)
summary(dun)
```

```
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Dunnett Contrasts
## 
## 
## Fit: aov(formula = VR ~ TRAT, data = dados_iva3)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value Pr(>|t|)
## B - A == 0    4.000      3.239   1.235    0.483
## C - A == 0    7.000      3.239   2.161    0.112
## D - A == 0   -4.000      3.239  -1.235    0.483
## (Adjusted p values reported -- single-step method)
```

``` r
confint(dun)
```

```
## 
##   Simultaneous Confidence Intervals
## 
## Multiple Comparisons of Means: Dunnett Contrasts
## 
## 
## Fit: aov(formula = VR ~ TRAT, data = dados_iva3)
## 
## Quantile = 2.5933
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##            Estimate lwr      upr     
## B - A == 0   4.0000  -4.4008  12.4008
## C - A == 0   7.0000  -1.4008  15.4008
## D - A == 0  -4.0000 -12.4008   4.4008
```

---

# Parte 8 — Fator quantitativo: regressão (solo)

---

## Dados (solo)

``` r
dados_solo <- tibble::tribble(
  ~TRAT, ~VR,
  1.31,2.61, 1.31,2.63, 1.31,2.65, 1.31,2.64, 1.31,2.62,
  1.43,2.57, 1.43,2.55, 1.43,2.59, 1.43,2.60, 1.43,2.56,
  1.55,2.50, 1.55,2.52, 1.55,2.48, 1.55,2.47, 1.55,2.46,
  1.67,2.42, 1.67,2.41, 1.67,2.39, 1.67,2.38, 1.67,2.40
) |>
  mutate(trat = factor(TRAT))
```

---

## ANOVA + diagnóstico (solo)

``` r
m_solo <- aov(VR ~ trat, data = dados_solo)
summary(m_solo)
```

```
##             Df  Sum Sq Mean Sq F value   Pr(>F)    
## trat         3 0.15274 0.05091   134.9 1.44e-11 ***
## Residuals   16 0.00604 0.00038                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

``` r
plot_diagnosticos(m_solo, "Diagnóstico — DIC Solo")
```

![](Análise-de-Variância-e-Planejamento-de-Experimentos_files/figure-html/reg-solo-diag-1.png)

``` r
shapiro.test(residuals(m_solo))
```

```
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(m_solo)
## W = 0.95503, p-value = 0.4499
```

``` r
bartlett.test(VR ~ trat, data = dados_solo)
```

```
## 
##  Bartlett test of homogeneity of variances
## 
## data:  VR by trat
## Bartlett's K-squared = 0.95836, df = 3, p-value = 0.8113
```

---

## Regressão linear + visual

``` r
m_lin  <- lm(VR ~ TRAT,             data = dados_solo)
m_quad <- lm(VR ~ TRAT + I(TRAT^2), data = dados_solo)

summary(m_lin)
```

```
## 
## Call:
## lm(formula = VR ~ TRAT, data = dados_solo)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.02920 -0.01415 -0.00250  0.01165  0.03860 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.48852    0.05074   68.75  < 2e-16 ***
## TRAT        -0.64833    0.03392  -19.12  2.1e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02035 on 18 degrees of freedom
## Multiple R-squared:  0.9531, Adjusted R-squared:  0.9504 
## F-statistic: 365.4 on 1 and 18 DF,  p-value: 2.1e-13
```

``` r
summary(m_quad)
```

```
## 
## Call:
## lm(formula = VR ~ TRAT + I(TRAT^2), data = dados_solo)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.03110 -0.01335 -0.00030  0.01335  0.03110 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   2.3416     0.6615   3.540  0.00252 **
## TRAT          0.9037     0.8934   1.012  0.32594   
## I(TRAT^2)    -0.5208     0.2996  -1.738  0.10023   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01929 on 17 degrees of freedom
## Multiple R-squared:  0.9601, Adjusted R-squared:  0.9554 
## F-statistic: 204.7 on 2 and 17 DF,  p-value: 1.273e-12
```

``` r
ggplot(dados_solo, aes(x = TRAT, y = VR)) +
  geom_point(size = 2, alpha = 0.8) +
  geom_smooth(method = "lm", formula = y ~ x,
              se = TRUE, linewidth = 1.2, colour = "#1B6CA8") +
  labs(title = "Regressão linear: densidade do solo vs MSPA",
       x = "Densidade (Mg/m³)", y = "MSPA (g)") +
  theme_minimal(base_size = 16)
```

![](Análise-de-Variância-e-Planejamento-de-Experimentos_files/figure-html/reg-solo-plot-1.png)

---

## Regressão com teste de desvios (ExpDes.pt)

``` r
if (requireNamespace("ExpDes.pt", quietly = TRUE)) {
  ExpDes.pt::dic(
    trat  = dados_solo$TRAT,
    resp  = dados_solo$VR,
    quali = FALSE,
    sigT  = 0.05,
    sigF  = 0.05
  )
} else {
  message("Instale o pacote: install.packages('ExpDes.pt')")
}
```

```
## ------------------------------------------------------------------------
## Quadro da analise de variancia
## ------------------------------------------------------------------------
##            GL      SQ       QM     Fc      Pr>Fc
## Tratamento  3 0.15274 0.050912 134.87 1.4392e-11
## Residuo    16 0.00604 0.000377                  
## Total      19 0.15877                           
## ------------------------------------------------------------------------
## CV = 0.77 %
## 
## ------------------------------------------------------------------------
## Teste de normalidade dos residuos ( Shapiro-Wilk ) 
## Valor-p:  0.4499191 
## De acordo com o teste de Shapiro-Wilk a 5% de significancia, os residuos podem ser considerados normais.
## ------------------------------------------------------------------------
## 
## ------------------------------------------------------------------------
## Teste de homogeneidade de variancia 
## valor-p:  0.8113263 
## De acordo com o teste de bartlett a 5% de significancia, as variancias podem ser consideradas homogeneas.
## ------------------------------------------------------------------------
## 
## Ajuste de modelos polinomiais de regressao
## ------------------------------------------------------------------------
## 
## Modelo Linear
## ==========================================
##    Estimativa Erro.padrao    tc    valor.p
## ------------------------------------------
## b0   3.4885     0.0484    72.0102     0   
## b1  -0.6483     0.0324    -20.0212    0   
## ------------------------------------------
## 
## R2 do modelo linear
## --------
## 0.990742
## --------
## 
## Analise de variancia do modelo linear
## ====================================================
##                      GL   SQ     QM     Fc   valor.p
## ----------------------------------------------------
## Efeito linear        1  0.1513 0.1513 400.85    0   
## Desvios de Regressao 2  0.0014 0.0007  1.87  0.18586
## Residuos             16 0.0060 0.0004               
## ----------------------------------------------------
## ------------------------------------------------------------------------
## 
## Modelo quadratico
## =========================================
##    Estimativa Erro.padrao   tc    valor.p
## -----------------------------------------
## b0   2.3416     0.6662    3.5151  0.0029 
## b1   0.9038     0.8997    1.0046  0.3301 
## b2  -0.5208     0.3017    -1.7263 0.1036 
## -----------------------------------------
## 
## R2 do modelo quadratico
## --------
## 0.998108
## --------
## 
## Analise de variancia do modelo quadratico
## ====================================================
##                      GL   SQ     QM     Fc   valor.p
## ----------------------------------------------------
## Efeito linear        1  0.1513 0.1513 400.85    0   
## Efeito quadratico    1  0.0011 0.0011  2.98  0.10355
## Desvios de Regressao 1  0.0003 0.0003  0.77  0.39454
## Residuos             16 0.0060 0.0004               
## ----------------------------------------------------
## ------------------------------------------------------------------------
## 
## Modelo cubico
## =========================================
##    Estimativa Erro.padrao   tc    valor.p
## -----------------------------------------
## b0  -8.3620     12.2513   -0.6825 0.5047 
## b1  22.6482     24.8681   0.9107  0.3760 
## b2  -15.1794    16.7560   -0.9059 0.3784 
## b3   3.2793     3.7479    0.8750  0.3945 
## -----------------------------------------
## 
## R2 do modelo cubico
## -
## 1
## -
## 
## Analise de variancia do modelo cubico
## ====================================================
##                      GL   SQ     QM     Fc   valor.p
## ----------------------------------------------------
## Efeito linear        1  0.1513 0.1513 400.85    0   
## Efeito quadratico    1  0.0011 0.0011  2.98  0.10355
## Efeito cubico        1  0.0003 0.0003  0.77  0.39454
## Desvios de Regressao 0    0      0      0       1   
## Residuos             16 0.0060 0.0004               
## ----------------------------------------------------
## ------------------------------------------------------------------------
```

---

# Encerramento

---

## Síntese

.pull-left[
**Boa análise começa no delineamento**
- randomização → independência
- replicação → estima o erro
- blocagem → reduz `$QM_E$`

**ANOVA é o teste global**
- qualitativo → médias / contrastes
- quantitativo → regressão + desvios
]

| Situação | Caminho |
|----------|---------|
| Fator qualitativo | Tukey / SNK / Dunnett |
| Fator quantitativo | Regressão → dose ótima |
| Gradiente no campo | DBC (bloco) |
| Dois gradientes | Quadrado Latino |
| F não significativo | Encerre aqui |
]

---

## Entregáveis PBL

.box-pbl[
**Caso Saúde (DIC vs DBC)**  
Justifique a escolha do delineamento e conclua sobre os protocolos.
]

.box-pbl[
**Caso Agronomia (DBC)**  
Aplique o Tukey e recomende o melhor tratamento com base nos grupos.
]

.box-pbl[
**Mini-relatório (1–2 p.)**  
Delineamento → ANOVA → Diagnóstico → Pós-ANOVA → Conclusão aplicada.
]

---

# Obrigado!

### Perguntas?