---
title: "Análise de Variância em blocos"
author: "Prof. João Tolentino"
date: "2025-05-14"
output:
  html_document:
    code_folding: hide
params:
  arquivo: NULL
  grupo: NULL
  bloco: NULL
  resposta: NULL
---

```{r, include=FALSE}
# opções
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(warning = FALSE)
``` 

```{r, include=FALSE}
# Carregar pacotes
if (!require(pacman)) install.packages("pacman")
pacman::p_load(
  "tidyverse", "readxl", "moments", "car", "gt",
  "emmeans", "multcomp", "multcompView"
)
```

# Análise da variável `r params$resposta`

### Importar o dados do arquivo
```{r}
extensao <- tolower(gsub(".*\\.", "", params$arquivo))
if (extensao == "xlsx") {
  dados <- readxl::read_excel(
    params$arquivo,
  )
} else if (extensao == "csv") {
  dados <- readr::read_csv(params$arquivo)
} else {
  stop("O arquivo não parece ser nem .xlsx nem .csv")
}
```

# Análise Exploratória

## Estatísticas descritivas por grupo

```{r}
dados |>
  group_by(.data[[params$grupo]]) |>
  summarise(
    n = n(),
    media = mean(.data[[params$resposta]]),
    desvpad = sd(.data[[params$resposta]]),
    var = var(.data[[params$resposta]]),
    .groups = "drop"
  ) |>
  gt(locale = "pt") |>
  fmt_integer(columns = n) |>
  fmt_number(columns = c("media", "desvpad", "var")) |>
  cols_align(columns = -1, "center")
```

## Boxplot comparativo entre os grupos

```{r}
dados |>
  ggplot(aes(
    x = fct_reorder(.data[[params$grupo]], .data[[params$resposta]], .fun = mean),
    y = .data[[params$resposta]]
  )) +
  geom_boxplot() +
  stat_summary(fun = "mean", shape = 23) +
  theme_classic()
```


## Estatísticas descritivas por bloco

```{r}
dados |>
  group_by(.data[[params$bloco]]) |>
  summarise(
    n = n(),
    media = mean(.data[[params$resposta]]),
    desvpad = sd(.data[[params$resposta]]),
    var = var(.data[[params$resposta]]),
    .groups = "drop"
  ) |>
  gt(locale = "pt") |>
  fmt_integer(columns = n) |>
  fmt_number(columns = c("media", "desvpad", "var")) |>
  cols_align(columns = -1, "center")
```

## Boxplot comparativo entre os blocos

```{r}
dados |>
  ggplot(aes(
    x = fct_reorder(.data[[params$bloco]], .data[[params$resposta]], .fun = mean),
    y = .data[[params$resposta]]
  )) +
  geom_boxplot() +
  stat_summary(fun = "mean", shape = 23) +
  theme_classic()
```



# Análise de variância
```{r}
formula_aov <- as.formula(paste(params$resposta, "~", params$grupo,"+", params$bloco))
mod1 <- lm(formula_aov, data = dados)
anova(mod1)
```

# Pressupostos

## Normalidade dos resíduos

### Gráfico Q-Q

```{r}
# Criar um data frame com os resíduos
residuos_df <- data.frame(residuos = mod1 |> residuals())

# Gráfico Q-Q com ggplot2
ggplot(residuos_df, aes(sample = residuos)) +
  stat_qq(color = "blue") + # Pontos dos quantis
  stat_qq_line(color = "red", linetype = "dashed") + # Linha de referência
  labs(
    title = "Gráfico Q-Q dos Resíduos",
    x = "Quantis Teóricos",
    y = "Quantis dos Resíduos"
  ) +
  theme_classic()
```

### Histograma dos resíduos com Índices descritivos
```{r}
## Calcular curtose, assimetria e teste de Shapiro-Wilk
curtose <- residuos_df |>
  moments::kurtosis()
assimetria <- residuos_df |>
  moments::skewness()
teste_shapiro <- residuos_df$residuos |>
  shapiro.test()


# Histograma com ggplot2

nbins <- nclass.Sturges(residuos_df$residuos) # Número de bins calculado pela regra de Sturges

ggplot(residuos_df, aes(x = residuos)) +
  geom_histogram(aes(y = after_stat(density)), bins = nbins, fill = "lightblue", color = "black") +
  geom_density(color = "red", linewidth = 1) +
  labs(
    title = "Histograma dos Resíduos",
    x = "Resíduos",
    y = "Densidade"
  ) +
  # Adicionar anotações de curtose, assimetria e p-valor do teste de Shapiro-Wilk
  annotate("text",
    x = Inf, y = Inf, hjust = 1.1, vjust = 1.5,
    label = paste("Assimetria:", round(assimetria, 2)),
    size = 4, color = "black"
  ) +
  annotate("text",
    x = Inf, y = Inf, hjust = 1.1, vjust = 3,
    label = paste("Curtose:", round(curtose, 2)),
    size = 4, color = "black"
  ) +
  annotate("text",
    x = Inf, y = Inf, hjust = 1.1, vjust = 4.5,
    label = paste("Shapiro-Wilk p-valor:", round(teste_shapiro$p.value, 4)),
    size = 4, color = "black"
  ) +
  theme_classic()
```




## Homogeneidade das variâncias
### Resíduos vs. ajustados

```{r}
## resíduos vs. ajustados
plot(mod1, which = 1)
```

### Boxplot dos resíduos por grupo

```{r}
## boxplot dos resíduos por grupo
dados |>
  ggplot(aes(x = .data[[params$grupo]], y = residuals(mod1))) +
  geom_boxplot() +
  theme_classic()
```

### Razão maior/menor desvio-padrão

```{r}
## razão maior/menor desvio-padrão
dados |>
  group_by(.data[[params$grupo]]) |>
  summarise(desvpad = sd(.data[[params$resposta]]), .groups = "drop") |>
  mutate(razao = max(desvpad) / desvpad) |>
  arrange(desc(razao)) |>
  gt(locale = "pt") |>
  fmt_number(columns = -1)
```

### Teste de Levene


```{r}
## teste de Levene

formula_levene <- as.formula(paste("residuals(mod1) ~", params$grupo))
car::leveneTest(formula_levene, data = dados)
```


## Maiores desvios dos resíduos

```{r}
## maiores desvios dos resíduos
residuals(mod1) %>%
  abs() %>%
  sort(dec = TRUE) %>%
  head()
```


