---
# System prepended metadata

title: CAA068 - Métodos Estatísticos Aplicados a Ciência Animal

---

# CAA068 - Métodos Estatísticos Aplicados a Ciência Animal
Segunda-feira - Sala de Multimídias do CLCE 8:20 - 12h20

Link para grupo de [Whatsapp](https://chat.whatsapp.com/HYCwCtGKKA8BSSbDAobunA?mode=ac_t)

Outros materiais para apoio
https://agronomiar.github.io/livroagro/an%C3%A1lise-conjunta.html
https://lec.pro.br/download/ivan/public/courses/caa110/main/
https://lec.pro.br/download/ivan/public/courses/caa068/main/

### Calendário 2025.2

| Data | Tema 2 | 
| -------- | -------- |
| 18/08/2025 |Apresentção da disciplina
| 25/08/2025 | Nivelamento em estatistica 
|01/09/2025 |Experimentação animal
|08/09/2025 |Princípios básicos de R e R-Studio
|15/09/2025|Análise Estatística Descritiva 
|22/09/2025 |Delineamentos experimentais
|29/09/2025 |Análise de variância (comparação de médias)
|06/10/2025| Relacionamento linear de variáveis - Análise de correlação 
|13/10/2025|Redução do número de variáveis - Análise de correspondência (ANACOR)
|20/10/2025| Análise de componentes principais (PCA)
|27/10/2025| Agrupamentos – Análise de cluster
|03/11/2025|Modelos de regressão (Machine Learning)
|10/11/2025|Atividade Avaliativa II
17/11/2025|Atividade Avaliativa III
|24/11/2025| Atividade Avaliativa IV

## Principais algoritimos usados em pesquisa biológica
![](https://hackmd.io/_uploads/r1MonDz63.png)

## qual teste estatístico usar
![](https://hackmd.io/_uploads/SkSjlKGph.jpg)

# Nivelamento em estatistica e instalação do R e Rstudio
[Nivelamento em estatistica (pdf)](https://docs.google.com/document/d/1rVW14UQtZyM4QpBD18Fuy0JFHc3-MNYu6bP-ck0VFlg/edit?usp=drive_link)

Links para download

[Software R](https://cran-r.c3sl.ufpr.br/) 

Escolher a opção do sistema operacional do computador
![](https://hackmd.io/_uploads/rJkRmVYa2.png)

Ir em instalar pela primeira vez
![](https://hackmd.io/_uploads/BkFkEEF63.png)

Escolher a versão mais atual
![](https://hackmd.io/_uploads/HJIz44Ka2.png)

Executar o instalador

Para instalar o Rstudio, vá para o site do [Rstudio](https://posit.co/download/rstudio-desktop/)

Escolher o sistema operacional e executar o instlador
![](https://hackmd.io/_uploads/S1HtNEtTn.png)



# Aula Experimentação animal



# Aula Princípios básicos de R e R-Studio
Os passos contidos nesse tutorial foram retirados da seguinte página com algumas modificações:
https://hackmd.io/@lamoroso92/minicurso23-r

## A interface do R
Constituida por algumas partes como:
1. **Script**: Local onde escrevemos o código (perene);
2. **Console**: Local onde escrevemos o código (temporário) e observamos os resultaydos;
3. **Ambiente/Histórico**: Variaveis/Objetos salvos; Histórico de códigos rodados;
4. **Arquivos/Plots/Pacotes/Ajuda**
![](https://i.imgur.com/8nxphp4.png)

## Teclas de Atalho

Atalho	                       Função
CTRL + ENTER        	Executa linha atual ou seleção
ALT + -               Adiciona o operador <-
CTRL + SHIFT + R    	Adiciona uma nova seção de código
CTRL + SHIFT + M    	Adiciona o operador %>%
CTRL + L            	Limpa o terminal do R

## Operadores 
**Operadores matemáticos:**
Permite realizar operações básicas de forma direta. Como uma calculadora

```{r}
# Operadores matemáticos:

2 + 3
2 * 3
2 - 3
2 / 3
2 ^ 3

```

**Operadores lógicos:**
Permite comparar valores e expressões, retorna valores boleanos: Verdadeiro (TRUE; T) ou Falso (FALSE; F).

```{r}
# Operadores lógicos:

2 == 3    # Igualdade
2 != 3    # Desigualdade
2 <= 3    # Menor (<); Menor ou igual (<=)
2 >= 3    # Maior (>); Maior ou igual (>=)

2 < 3 & 3 > 5    # & = E; Só TRUE quando ambos são TRUE
2 < 3 | 3 > 5    # | = OU; É TRUE quando ao menos um é TRUE

```

**Operadores de atribuição:**
Permite atribuir valores a um objeto. São válidos os símbolos ‘<-’ e ‘=’.

```{r}
# Operadores de atribuição:

a <- 3   # Cria um objeto 'a' de valor 3
a        # Visualiza o objeto 'a'

b = 2     
b

a + b    # Operações válidas para o valor interno serão executadas normalmente.
a - b
a == b

```

## Vetores
Conjunto de elementos de um mesmo tipo (numérico, textual, lógicos…).

Criação de vetores:

```{r}
# Criação de vetores:

# Criando vetores com a função 'c()' (Concatenar).

v1 <- c(2, 6, 30, 193)    # Vetor numérico (Numeric)
v1

v2 <- c("a", "z", "lucas", "e. coli")    # Vetor de texto (Character)
v2

v3 <- c(TRUE, FALSE, TRUE)    # Vetor lógico (Logical)
v3

# Criando vetores de outras formas:

v <- 1:5    # Sequência de 1 a 5, o mesmo que c(1, 2, 3, 4, 5)
v

v <- rep(1, 7)    # Função rep(); Repetirá 1, 7 vezes
v

v <- rep("a", 7)    # Também funciona com texto
v

v <- seq(10, 20, 2)    # Função seq(); Sequência de 10 a 20, pulando de 2 em 2.
v
```

**Operação com vetores**

```{r}
# Operação com vetores

(vx <- 1:5)
(vy <- rep(1, 5))

# Operações podem ser feita entre os elementos, de forma pareada, quanto possuirem o mesmo tamanho:

length(vx)    # Função length() mede o tamanho de um objeto
length(vy)

vx + vy
vx + rep(1, 6)    # Erro
vx > vy   # Checa cada elemento de 'vx' contra cada elemento de 'vy' (na mesma posição)


# Ou para cada elemento:

vx * 2    # Multiplica cada elemento por 2
vx > 2    # Todos contra 2 

```

**Elementos de vetores**
```{r}
# Elementos de vetores

vl <- letters    # Objeto especial, contém alfabeto minúsculo
vl

vl[20]    # 20º elemento
vl[c(15, 9)]    # 15º e 9º elemento
vl[1:5]    # Série de elementos do 1º ao 5º

vn <- 1:10
vn

vn[c(T, T, F, F, T, T, F, F, T, T)]    # Seleção lógica; Exibe T
vn[vn > 5]    # Maiores que 5
vn[(vn %% 2) == 0]    # Apenas pares
vn[(vn %% 2) != 0]    # Apenas ímpares

```

## Matrizes e Data Frames
**Matriz:** Vetores com duas dimensões (linhas/rows e colunas/columns). Mesmo tipo!
```{r}
# Matrizes 

# Criação de matrizes:

m <- 1:20    # Vetor
m

dim(m)    # Função que descreve n. de dimensões

dim(m) <- c(5, 4)    # Alterando n. de dimensões
m

dim(m) <- c(10, 2)    # Alterando n. de dimensões
m

t(m)    # Invertendo (transpondo) a matriz
m       # Note que as mudanças não são perenes...

m <- t(m)   # ...A não ser que sejam salvas
m

# Seleção de elementos:

m[1, ] # 1ª linha
m[ ,5] # 5ª coluna
m[2,3] # Elemento da 2ª linha e 3ª coluna
m
```

**Data Frames:** Tabelas contendo vetores ou fatores de mesmo tamanho. Pode possuir tipos diferentes. Linhas/Rows também podem ser chamadas de Casos/Cases e Colunas/Columns de Váriaveis/Variable.
```{r}
# Data Frames

# Criação de DFs

df1 <- data.frame(
    ID = 1:5,
    Nome = c("Kuririn", "Goku", "Cell", "Piccolo", "Bulma"),
    Sexo = as.factor(c("M", "M", "M", "M", "F")),
    Idade = c(42, 46, 6, 26, 50),
    Especie = as.factor(c("Humano", "Sayajin", "Android", "Namekusei", "Humano"))
)

df1

# Seleção de elementos:

df1[1, ] # 1ª linha, todas as colunas
df1[ ,2] # 2ª coluna, todas as linhas

df1$Nome    # Coluna "Nome"
df1["Nome"] # Coluna "Nome"

df1[c(1,3), ] # Linhas 1 e 3, todas as colunas
df1[1:2, 2:5] # Linhas 1 e 2, colunas de 2 a 5
df1[df1$Idade > 30, c("Nome", "Especie")] # Das linhas cuja idades seja maior que 30, mostrar colunas "Nome" e "Especie"

# Adicionando colunas

df1$Poder <- c(150, 9000, 5500, 300, 10)
df1
```

## Controle de fluxo
Processos iterativos (em loop) que “quebram” o fluxo do código, realizando tarefas até que uma determinada condição seja satisfeita.

for
Para cada elemento, faça algo.
```{r}
# for:

# Para cada elemento de 1 a 10
for (i in 1:10) {
    
    # Exiba o resultado do elemento * 2
    print(i * 2)
    
}

```
**if/else**
Se uma condição for verdadeira, faça algo; Se não, faça outra coisa.
```{r}
# if/else:

x <- 4

# Se módulo de x/2 for igual a 0
if (x %% 2 == 0) {
    
    # Exiba "par"
    print("par")

# Se não
} else {
    
    # Exiba "ímpar"
    print("ímpar")
    
}

```

**while**
Enquanto uma condição se manter verdadeira, faça algo.
```{r}
# while:

x <- 1

# Enquanto 'x' for menor que 10, faça:
while (x <= 10) {
    
    # Exiba:
    print(paste("'x' agora é", x))
    
    # Some 1 ao atual valor de 'x'
    x <- x + 1
    
}
```

**Funções**
Aplicações para realizar ações específicas. Quase tudo é uma função. Das mais simples às mais complexas.

Algumas funções já utilizadas: *c(), length(), rep(), seq(), t(), data.frame()*…

Funções são pertencentes a pacotes nativos (built-in), carregados (instalados e convocados) ou criadas.

Funções podem receber *argumentos* limitados, ilimitados, ou nem serem necessários. Alguns são opcionais, alguns são obrigatórios.

Argumentos pode ou não ser nomeados. Se não forem, devem respeitar a ordem estabelecida pela função.

Para checar os requisitos, usabilidade e exemplos do uso de uma função, utilize um ? antes de seu nome. (Ex. ?length)

**Built-in**
As funções mais básicas são nativas. Pertencem a pacotes iniciados em conjunto com o R.
```{r}
# Funções built-in:

# Funções sem necessidade de argumentos:

# ls(): Lista os objetos
ls()

# getwd(): Informa o local de trabalho
getwd()

# Argumentos nomeados vs. não nomeados:
?round()
round(x = 1.7468357)
round(x = 1.7468357, digits = 3)
round(digits = 3 , x = 1.7468357)

round(1.7468357, 3)
round(3, 1.7468357)

# Funções específicas para determinados tipos de valores:
sum(23, 12, 34, 45, 12, 59)
sum("a", "b", "a", "d", "f") # erro

```


**Funções criadas**
É possivel criar uma função nova com a função: function().
```{r}
# Funções criadas:

critico <- function(titulo, nota){
    
    if (nota <= 4 ){
        paste("O filme", titulo, "é RUIM! Guarde seu dinheiro.")
    } else if (nota > 4 & nota < 8){
        paste("O filme", titulo, "é MÉDIÃO! Só veja se não tiver algo melhor pra fazer.")
    } else {
        paste("O filme", titulo, "é ÓTIMO! Pode ver sem medo.")
    }
    
}

critico(titulo = "Morbius", nota = 3.5)
critico("Avatar", 5.5)
critico(nota = 9.5, titulo = "O Senhor dos Anéis")

```

## Pacotes
Pacotes são conjuntos de funções que ampliam as capacidades do R. Esses pacotes devem ser instalados a parte e não são carregados automaticamente, requerendo uma convocação em seu primeiro uso.

O R não possui a capacidade nativa de lidar com arquivos no formato excel. O pacote xlsx fornece funções para abrir e salvar arquivos nesse formato.
[Podekex](https://docs.google.com/spreadsheets/d/10RFtI26jeuOcPBvmTysHqWMnJX6Uz7_o/edit?usp=sharing&ouid=116739411541654462827&rtpof=true&sd=true)

Library xlsx precisa do [java](https://www.java.com/en/download/)

```{r}
# Funções de pacotes:

# Instalação de pacotes - Execução única!
install.packages("xlsx")
install.packages("readxl")


# Convocação do pacote - Uma vez por sessão
library("xlsx") # require("xlsx") também funciona
library("readxl")

# Funções de um pacote:
# Digitar "pacote::" + tecla 'tab', ou:
ls("package:xlsx")

# Buscando ajuda:
?read.xlsx

# Abrindo tabela do excel:
pokedex <- read.xlsx(file = "/caminho/para/arquivo/pokedex.xlsx", sheetIndex = 1)
pokedex <- read_excel("D:/pokedex.xlsx")


# Visualizando linhas iniciais, função :
head(pokedex)

# Visualizadno tabela inteira:
View(pokedex)

```

## Tidyverse
Coleção de pacotes do R;
Compartilham estilo e síntaxe; Grande parte criado por Hadley Wickham;
Voltado à processos de data wrangling:
Obtenção, transformação, limpeza, agregação e visualização dos dados
Além de pacotes como dplyr, tidyr e ggplot2, o tidyverse adiciona o operador pipe (%>%), que permite o encadeamento de funções de forma mais intuitiva:

```{r}
# Instalação do 'tidyverse' já foi realizada préviamente (pesada!).
# install.packages("tidyverse")

# Convocação
install.packages("tidyr")
install.packages("dplyr")
library("tidyr") # pacote tidyr
library("dplyr") # pacote dplyr
library("tidyverse")

# Função única:
sum(13, 18, 36, 16)

sqrt(83)

round(9.110434, digits = 2)

# Ou

a <- sum(13, 18, 36, 16)

b <- sqrt(a)

round(b, digits = 2)

# Método de encadear funções no R Base (aninhamento):
round(x = sqrt(sum(13, 18, 36, 16)), digits = 2)

# Método com o pipe "%>%"
sum(13, 18, 36, 16) %>% sqrt() %>% round(digits = 2)

```
## Básicos da exploração de dados em R

**Exploração**

Estrutura dos dados

```{r}
# Recarregando a tabela:
pokedex <- read_excel("D:/pokedex.xlsx")

# Visualização completa da tabela (aba):
View(pokedex)

# Visualização das primeiras linhas:
head(pokedex)

# Tamanho da tabela:
ncol(pokedex)     # N. de colunas
nrow(pokedex)     # N. de linhas

dim(pokedex) # N. linhas e N. colunas 

# Resumo da tabela (tipos):
str(pokedex)

```
Variáveis, tipos e valores únicos

```{r}
# Checando os nomes:
names(pokedex)

# Tipo de variáveis:
class(pokedex$Name)
class(pokedex$Height_m)

# Valores únicos:
unique(pokedex$Type1)

```
Dados faltantes
```{r}
# Dados faltantes:
# Ausência de dados pode ser um problema. No R, células vazias recebem 'NA' (not available).
pokedex[ complete.cases(pokedex), ]
pokedex[!complete.cases(pokedex), ]
```

**Limpeza**
Seleção de variáveis
```{r}
# Seleção de variáveis. Função 'select':
names(pokedex)

pokedex %>% 
  select(1:2, 4) # por 'posição'

pokedex %>% 
  select(Name:Type1, Attack) # por 'nomes'

pokedex %>% 
  select(-Type1, -Speed, -`Sp. Atk`, - Weight_kg, -Height_m) # seleção negativa

pokedex %>% 
  select(Name, starts_with("Type")) # por padrão

pokedex2 <- pokedex %>% 
  select(Name:Defense, `Sp. Atk`:Height_m) %>% 
  select(-`Sp. Atk`, -Speed)

pokedex2

```

Mudando ordem de variáveis
```{r}
# Mudando ordem de variáveis. Função 'select':
pokedex2 %>% 
    select(Name, everything())
```

Mudando nomes de variáveis
```{r}
# Mudando nome de variáveis. Função 'rename':
pokedex2 <- pokedex2 %>%
    # Sintaxe: rename(nome_novo = nome_velho)
    rename(Type = Type1)
```

Mudando o tipo de variáveis
```{r}
# Mudando o tipo de variáveis:
class(pokedex2$Type)

pokedex2$Type <- as.factor(pokedex2$Type)

class(pokedex2$Type)

levels(pokedex2$Type)
```

Filtrando linhas

```{r}
# Filtrando linhas:
pokedex2 %>%
    select(Name, Type, Attack) %>%
    # Filtando por uma condição
    filter(Type == "Fire")

pokedex2 %>%
    select(Name, Type, Attack) %>%
    # Filtando se uma condição OU outra for verdadeira
    filter(Type == "Fire" | Type == "Water")

pokedex2 %>%
    select(Name, Type, Attack) %>%
    # Filtando se uma condição E outra forem verdadeiras
    filter(Type == "Fire" & Attack >= 100)

pokedex2 %>%
    select(Name, Type, Attack) %>%
    # Soma dos filtros anteriores
    filter((Type == "Fire" | Type == "Water") & Attack >= 100)
    
```

Lidando com dados perdidos (“NAs”)
```{r}

# Transformando a coluna em valores numericos
pokedex$Weight_lb <- as.numeric(pokedex$Weight_lb)

# Lidando com "NAs"
mean(pokedex$Weight_lb ) # Erro
mean(pokedex$Weight_lb, na.rm = TRUE)

pokedex2 <- pokedex %>%
    drop_na()

mean(pokedex2$Weight_lb)
```
**Manipulação**

Criando ou alterando variável

```{r}
# Criando ou alterando variáveis. Função 'mutate':
pokedex2 <- pokedex2 %>%
    mutate(Height_m = Height_cm / 100) %>%
    mutate(Weight_kg = Weight_lb / 2.205) %>%
    select(-Height_cm, -Weight_lb)

View(pokedex2)

```

Mudança condicional
```{r}
# Mudança condicional. Funtção 'mutate' + 'if_else':
pokedex2 <- pokedex2 %>%
    mutate(Power = if_else(Attack >= 120 | Defense >= 60, "High", "Low"))

View(pokedex2)

```
**Descrição**

Amplitude e Dispersão
```{r}
# Amplitude: Mínimas e máximas:
min(pokedex2$Height_m)
max(pokedex2$Height_m)
range(pokedex2$Height_m)

# Dispersão
IQR(pokedex2$Height_m)
boxplot(pokedex2$Height_m, horizontal = T)

```
Centralidade
```{r}
# Centalidade
mean(pokedex2$Attack)
median(pokedex2$Attack)

```
Variância e Desvio-Padrão
```{r}
# Variância e Desvio-Padrão
var(pokedex2$Attack)
var(pokedex2$Height_m)

sd(pokedex2$Attack)
sd(pokedex2$Height_m)

## grade de grafico 1x2
par(mfrow=c(1,2))
boxplot(pokedex2$Attack, horizontal = F)
boxplot(pokedex2$Height_m, horizontal = F)


## grade de grafico 1x1 
par(mfrow=c(1,1))
boxplot(pokedex2$Attack, horizontal = F)
boxplot(pokedex2$Height_m, horizontal = F)

```
Resumo dos dados
```{r}
# Resumos
summary(pokedex2$Attack)
summary(pokedex2$Height_m)

pokedex2 %>% 
    select(Attack, Defense, Height_m, Weight_kg) %>%
    summary()

```
Criação de tabelas de resumos
```{r}
# Criação de tabelas de resumos
table(pokedex2$Type)

tipo_freq <- pokedex2 %>%
    select(Type) %>%
    table()

View(tipo_freq)

tipo_atk <- pokedex2 %>%
    # Remove colunas não-agrupaveis e não-numéricas
    select(-Name) %>%
    # Agrupa pelo fator em comum
    group_by(Type) %>%
    # Cria as medidas as quais se deseja avaliar
    summarise(Min = min(Attack),
              Max = max(Attack),
              Mean = mean(Attack),
              Var = var(Attack),
              Amplitude = max(Attack)-min(Attack),
              SD = sd(Attack)) %>%
    # Ordena de forma decrescente ('-')
    arrange(-Mean)

tipo_df <- pokedex2 %>%
    select(-Name) %>%
    group_by(Type) %>%
    summarise(Min = min(Defense ),
              Max = max(Defense ),
              Mean = mean(Defense ),
               Var = var(Attack),
              Amplitude = max(Attack)-min(Attack),
              SD = sd(Defense )) %>%
    arrange(-Mean)
    
    
    
#### criar um resumo para a planilha com os avlores numericos 
tipo_resumo <- pokedex2 %>% 
  select(-Name) %>% 
  group_by(Type1) %>% 
  summarise(
    across(
      where(is.numeric),
      list(
        Min  = ~min(.x, na.rm = TRUE),
        Max  = ~max(.x, na.rm = TRUE),
        Mean = ~mean(.x, na.rm = TRUE),
        Var  = ~var(.x, na.rm = TRUE),
        SD   = ~sd(.x, na.rm = TRUE)
      ),
      .names = "{.col}_{.fn}"
    )
  )
    

View(tipo_atk)
View(tipo_df)

```
## Visualização
ggplot2

Um dos melhores pacotes do R para visualização de dados;
Permite a construção modular de gráficos;
Possui uma “gramática” própria (“A Gramática dos Gráficos”). Básicamente, necessita de três informações: data, mapping e geometry
Inicialmente parece complexo, mas é compensado por sua versatilidade;
Cheat Sheet: Link

Gráficos de barra
```{r}
# Gráfico de barras
ggplot(data = pokedex2, mapping = aes(x = Type)) +
    geom_bar()

ggplot(data = pokedex2, mapping = aes(x = Type, fill = Type)) +
    geom_bar()

cores <- c("#a8b821", "#6e39fd", "#f7d432", "#c03121",
           "#f07f2f", "#6f5997", "#78c750", "#e0c167",
           "#98d8d8", "#a9a878", "#9f3fa0", "#fc5686",
           "#b89f38", "#6294e9","#d25252")

ggplot(data = pokedex2, mapping = aes(x = Type, fill = Type)) +
    geom_bar() +
    scale_fill_manual(values=cores)
    
    
ggplot(data = pokedex2, mapping = aes(x = Type1, fill = Type1)) +
  geom_bar() +
  labs(
    title = "Distribuição dos Tipos de Pokémon",
    x = "Tipo Primário",
    y = "Número de Pokémon",
    fill = "Legenda de Tipos"
  ) +
  theme_classic() +
  #scale_fill_manual(values=cores)+
  theme(
    plot.title = element_text(size = 14, face = "bold", hjust = 0.0), # título centralizado
    axis.title.x = element_text(size = 14, face = "bold"),
    axis.title.y = element_text(size = 14, face = "bold"),
    axis.text.x = element_text(size = 14, angle = 45, hjust = 1), # texto do eixo X inclinado
    axis.text.y = element_text(size = 14),
    legend.title = element_text(size = 14, face = "bold"),
    legend.text = element_text(size = 14)
  )
    
    

```

Histogramas
```{r}
# Histogramas
ggplot(data = pokedex2, mapping = aes(x = Height_m)) +
    geom_histogram()

ggplot(data = pokedex2, mapping = aes(x = Height_m)) +
    geom_histogram(colour = "black", fill = "white", bins = 100)

ggplot(data = pokedex2, mapping = aes(x = Height_m)) +
    geom_histogram(colour = "black", fill = "white", bins = 20)

```
Diagrama de caixa (Boxplot)
```{r}
# Diagrama de caixa (Boxplot)
ggplot(data = pokedex2, mapping = aes(y = Defense)) +
    geom_boxplot()

ggplot(data = pokedex2, mapping = aes(y = Attack, x = Type)) +
    geom_boxplot()

ggplot(data = pokedex2, mapping = aes(y = Speed, x = Type, fill = Type)) +
    geom_boxplot() +
    scale_fill_manual(values=cores) +
    theme_bw()
```

Gráfico de densidade
```{r}
# Gráfico de densidade
pokedex2 %>%
    filter(Type %in% c("Grass", "Fire")) %>%
    ggplot(aes(x = Attack, color = Type, fill = Type)) +
    geom_density(alpha = 0.8) +
    theme_bw() +
    scale_fill_manual(values=c("#c03121", "#fc5686")) +
    scale_color_manual(values=c("#c03121", "#fc5686")) 

pokedex2 %>%
    filter(Type %in% c("Bug", "Water")) %>%
    ggplot(aes(x = Defense, color = Type, fill = Type)) +
    geom_density(alpha = 0.8) +
    theme_bw() +
    scale_fill_manual(values=c("#c03121", "#fc5686")) +
    scale_color_manual(values=c("#c03121", "#fc5686")) 

```
Gráfico de dispersão
```{r}
# Gráfico de dispersão
pokedex2 %>%
    ggplot(aes(Speed, HP)) +
    geom_point(size = 3, alpha = 0.8) +
    theme_bw() +
    labs(title = "Velocidade vs. HP")

pokedex2 %>%
    ggplot(aes(Attack, Defense)) +
    geom_point(size = 3, alpha = 0.8) +
    theme_bw() +
    labs(title = "Defesa vs. Ataque")

pokedex2 %>%
    ggplot(aes(Attack, Defense, color = Type)) +
    geom_point(size = 3, alpha = 0.8) +
    theme_bw() +
    scale_color_manual(values=cores) +
    labs(title = "Defesa vs. Ataque")

```
Gráfico de tendência
```{r}
# Gráfico de tendência
pokedex2 %>%
    ggplot(aes(Attack, Defense)) +
    geom_point(size = 3, alpha = 0.8) +
    geom_smooth() +
    theme_bw() +
    labs(title = "Defesa vs. Ataque")

pokedex2 %>%
    ggplot(aes(Attack, Defense)) +
    geom_point(size = 3, alpha = 0.8) +
    geom_smooth(method = "lm") +
    theme_bw() +
    labs(title = "Defesa vs. Ataque")

```

# Análise Estatística Descritiva 
As etapas contidas nesta seção foram retiradas das seguintes páginas com modificações:
https://rpubs.com/paulofilho/751903
https://agronomiar.github.io/livroagro/estat%C3%ADstica-descritiva.html

```{r}

# instalando vários pacotes e checando se estão instalados
# Lista de pacotes
pacotes <- c("tidyverse", "DT", "knitr", "maps", "mapdata", "devtools")

# Identifica pacotes que não estão instalados
nao_instalados <- setdiff(pacotes, rownames(installed.packages()))

# Instala apenas os que faltam
if(length(nao_instalados) > 0){
  install.packages(nao_instalados, dependencies = TRUE)
}

# Carrega todos
sapply(pacotes, require, character.only = TRUE)

```

## Medidas de tendência central
As principais medidas usadas para averiguar a tendência central são: média, mediana e moda. 

### Média
Observando os dados

```{r}
# Banco de dados Iris do próprio R

head(iris) # primeiras 6 linhas do dataframe
datatable(head(iris, 20)) # tabela mais bonita e interativa
# colocando somente os valores numéricos
vars <- iris[, -5] # todos as colunas, exceto a última


```

A média aritmética simples é A medida de tendência central mais comumente usada para descrever resumidamente um conjunto de dados. É definida como a soma das observações dividida pelo número delas. 

A função tapply em R é usada para aplicar uma função específica a subconjuntos de dados de acordo com os níveis de um ou mais fatores. Ela é frequentemente usada para realizar operações agregadas em dados divididos por categorias ou fatores.


```{r}
# Algumas formas de calcular a média do dataframe

# usando a função tapply
tapply(X = iris$Sepal.Length, INDEX = list(iris$Species), FUN = mean)

```
A função aggregate em R é usada para realizar operações de agregação em dados, geralmente para resumir informações em grupos com base em fatores ou variáveis. Ela é semelhante à função tapply, mas oferece mais flexibilidade e é frequentemente usada quando você precisa de um controle mais detalhado sobre como os dados são agregados

```{r}
# Algumas formas de calcular a média do dataframe

# usando a função aggregate
aggregate(Sepal.Length ~ Species, data = iris, mean)

```
A função apply em R é usada para aplicar uma função a matrizes ou arrays ao longo de dimensões específicas


```{r}

# usando a função apply

apply(vars, 2, mean) # usa o objeto "var", na coluna [2] e aplica a média (mean)

```

### Mediana
A mediana procura caracterizar o centro da distribuição de frequências quando os valores são dispostos em ordem crescente ou decrescente de magnitude.
É o valor que divide o conjunto ordenado de valores em duas partes com igual número de elementos, ou seja, 50% das observações ficam acima da mediana e 50% ficam abaixo.

```{r}
# Algumas formas de calcular a mediana do dataframe

# função aggregate
aggregate(Sepal.Length ~ Species, data = iris, median)

# função apply
apply(vars, 2, mean) # usa o objeto "var", na coluna [2] e aplica a mediana (mean)

# função tapply
tapply(X = iris$Sepal.Length, INDEX = list(iris$Species), FUN = mean)

```
### Moda
Valor mais frequente

```{r}
moda <- sort(table(iris$Petal.Length), decreasing = TRUE)
moda[1]
```

### Máximo

Valor máximo observado 
(máximo = max(iris$Petal.Length))
```{r}
maximo = max(iris$Petal.Length)
maximo
```

### Mínimo
Valor mínimo observado 
```{r}
minimo = min(iris$Petal.Length)
minimo
```

## Medidas de dispersão 
As medidas de dispersão servem para indicar o quanto os dados se apresentam dispersos, ou afastados, em relação ao seu valor médio.

### Amplitude Total
A maneira mais simples de se medir a variabilidade de uma variável é através da “distância” entre o maior e o menor valor observado em um conjunto de dados.

```{r}
amplitude = max(iris$Petal.Length) - min(iris$Petal.Length)
amplitude

```
### Variância da amostra
A variância mede o quanto os dados estão dispersos em torno da média, ou seja, do valor esperado.

É a soma dos quadrados dos desvios, dividida pelo total de observações menos um.

```{r}
apply(vars, 2, var)
tapply(X = iris$Sepal.Length, INDEX = list(iris$Species), FUN = var)
aggregate(Sepal.Length ~ Species, data = iris, var)

```

### Desvio Padrão
O desvio padrão é uma medida que indica a dispersão dos dados dentro de uma amostra com relação à média.
É a raiz quadrada da variância

```{r}
apply(vars, 2, sd)
tapply(X = iris$Sepal.Length, INDEX = list(iris$Species), FUN = sd)
aggregate(Sepal.Length ~ Species, data = iris, sd)
```

### Coeficiente de Variação
O coeficiente de variação de Pearson é a razão entre o desvio padrão e a média. Em geral, o resultado é multiplicado por 100, para que o coeficiente de variação seja expresso em porcentagem.

```{r}
CV = sd(iris$Petal.Length) / mean(iris$Petal.Length)*100
CV
```

### gerando tabelas simplificadas

```{r}
by(iris, iris$Species, summary)

```

### Criando novas caracteristicas

Criar uma nova coluna "Size" na qual se o "Sepal.Length" for menor que a média será "small" e se for maior "big"

```{r}
iris$size <- ifelse(iris$Sepal.Length < median(iris$Sepal.Length),
                    "small", "big")
```

### Tabela de contingência 
Pode-se criar uma tabela de contingência baseado nos para saber a relação entre as espécies e o tamanho do septo das folhas
```{r}

table(iris$Species, iris$size)
```

### Correlação
A função cor em R é usada para calcular a correlação entre variáveis numéricas. A correlação é uma medida estatística que indica o grau de relacionamento linear entre duas variáveis. Ela ajuda a determinar se e como duas variáveis estão relacionadas entre si.

```{r}
cor(vars) # todos os valores 

cor(iris$Sepal.Length, iris$Sepal.Width) # valores especificos

```

### Visualização em gráficos

#### O gráfico em mosaico
permite observar a tabela de contingência de forma simplificada

```{r}
mosaicplot(table(iris$Species, iris$size),
           color = TRUE,
           xlab = "Species", 
           ylab = "Size" 
)
```

#### Grafico em barras 
Um barplot é uma ferramenta para visualizar a distribuição de uma variável qualitativa.

```{r}
barplot(table(iris$size)) # table() is mandatory

barplot(prop.table(table(dat$size)))

ggplot(iris) +
  aes(x = size) +
  geom_bar()+
  theme_classic()


```

#### Histograma 
Um histograma fornece uma ideia sobre a distribuição de uma variável quantitativa.

```{r}
hist(iris$Sepal.Length)
ggplot(iris) +
  aes(x = Sepal.Length) +
  geom_histogram()
  
  ```
  
#### Boxplot
Os boxplots são realmente úteis na estatística descritiva e muitas vezes são subutilizados (principalmente porque não são bem compreendidos pelo público). Um boxplot representa graficamente a distribuição de uma variável quantitativa, exibindo visualmente cinco resumos comuns de localização (mínimo, mediana, primeiro e terceiro quartis e máximo) e qualquer observação que tenha sido classificada como um valor atípico suspeito usando o critério do intervalo interquartil (IQR).

O critério do IQR significa que todas as observações acima de q0.75 + 1.5 * IQR ou abaixo de q0.25 - 1.5 * IQR (onde q0.25 e q0.75 correspondem ao primeiro e terceiro quartis, respectivamente) são consideradas como possíveis valores atípicos pelo R. O mínimo e o máximo no boxplot são representados sem esses valores atípicos suspeitos.

Ver todas essas informações no mesmo gráfico ajuda a ter uma boa primeira visão geral da dispersão e localização dos dados. Antes de traçar um boxplot dos nossos dados, veja abaixo um gráfico explicando as informações presentes em um boxplot.

![](https://hackmd.io/_uploads/BJjEyxlJ6.png)
fonte: https://statsandr.com/blog/descriptive-statistics-in-r/


```{r}
boxplot(iris$Sepal.Length)
boxplot(iris$Sepal.Length ~ iris$Species)
ggplot(iris) +
  aes(x = Species, y = Sepal.Length) +
  geom_boxplot()
  
```

#### Dotplot
Um dotplot é mais ou menos semelhante a um boxplot, exceto que:

- As observações são representadas como pontos.
- Ele não fornece facilmente informações sobre a mediana, primeiro quartil e terceiro quartil.




```{r}
dotplot(iris$Sepal.Length ~ iris$Species)

# vertical
ggplot(iris) +
  aes(x = Species, y = Sepal.Length) +
  geom_dotplot(binaxis = "y", stackdir = "center")

# horizontal
ggplot(iris) +
  aes(x = Sepal.Length, y = Species) +
  geom_dotplot(binaxis = "y", stackdir = "center")

```

### Scatterplot
Scatterplots permitem verificar se existe uma possível ligação entre duas variáveis quantitativas.

```{r}

plot(iris$Sepal.Length, iris$Petal.Length)

ggplot(iris) +
  aes(x = Sepal.Length, y = Petal.Length) +
  geom_point()

ggplot(iris) +
  aes(x = Sepal.Length, y = Petal.Length, colour = Species) +
  geom_point() +
  scale_color_hue()
```

### Lineplot
Gráficos de linhas, particularmente úteis em séries temporais ou em relação linear das variaveis, podem ser criados adicionando o argumento type = "l" na função plot():
```{r}
plot(dat$Sepal.Length,
  type = "l"
) # "l" for line
```




### Relação entre 2 ou mais variáveis
Se quisermos visualizar vários planos formados por pares de vetores numéricos, podemos gerar um painel dos gráficos de dispersão com a função pairs():
```{r}
pairs (iris[,1:4], 
       col = as.numeric(iris$Species), 
       oma = c(3,3,3,14))
par (xpd = TRUE)
legend ("bottomright",
        pch = c(1,1,1),
        col = unique(as.numeric(iris$Species)), 
        legend = levels(iris$Species))
```

# Delineamento experimentais + Análise de variância (comparação de médias)
[Aula]
[anova](https://drive.google.com/file/d/10g2oFlc-JmVS42znlNeiTibWBCv-EFGK/view?usp=sharing)
[bloco](https://drive.google.com/file/d/10hD9Dpoi0rkAL5YUivxG3ixyQRQsKDaA/view?usp=sharing)

## analise paramétrica
```{r}
library("dplyr")
library("ggpubr")
library("car")
library("agricolae")
library("DescTools")



dados <- read.csv("D:/anova.csv")
str(dados)

# Criando os vetores com os dados
Tratamento <- c("A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "C", "C", "C", "C", "C", "D", "D", "D", "D", "D")
Leite <- c(23, 21, 22, 20, 21, 25, 22, 23, 20, 19, 26, 27, 20, 24, 22, 27, 28, 25, 23, 24)
Bloco <- c("B1","B2","B3","B4","B5","B1","B2","B3","B4","B5","B1","B2","B3","B4","B5",
           "B1","B2","B3","B4","B5")
# Criando o dataframe
dados <- data.frame(Tratamento, Leite, Bloco)

# Visualizando o dataframe
print(dados)


dados
str(dados)
levels(dados$Tratamento)
# dados$Tratamento <- as.factor(dados$Tratamento)
# dados$Tratamento <- as.character(dados$Bloco)

levels(dados$Tratamento)


group_by(dados, Tratamento) %>%
  summarise(
    count = n(),
    mean = mean(Leite, na.rm = TRUE),
    sd = sd(Leite, na.rm = TRUE)
  )

ggboxplot(dados, x = "Tratamento", y = "Leite", 
          color = "Tratamento",
          order = c("A", "B", "C", "D"),
          ylab = "Leite", xlab = "Tratamento")

ggboxplot(dados, x = "Bloco", y = "Leite", 
          color = "Bloco",
          order = c("B1", "B2", "B3", "B4", "B5"),
          ylab = "Leite", xlab = "Bloco")

# One way anova DIC
anova <- aov(Leite ~ Tratamento, data = dados)
summary(anova)

# Two way anova DBC
anova2 <- aov(Leite ~ Tratamento+Bloco, data = dados)
summary(anova2)

#interação tratamento * bloco
anova3 <- aov(Leite ~ Tratamento*Bloco, data = dados)
summary(anova3)


#interacao com fatorial 
Tratamento_2 <- c("A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "C", "C", "C", "C", "C", "D", "D", "D", "D", "D",
                "A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "C", "C", "C", "C", "C", "D", "D", "D", "D", "D")
Leite_2 <- c(23, 21, 22, 20, 21, 25, 22, 23, 20, 19, 26, 27, 20, 24, 22, 27, 28, 25, 23, 24,
           23, 21, 22, 20, 21, 25, 22, 23, 20, 19, 26, 27, 20, 24, 22, 27, 28, 25, 23, 24)
Bloco_2 <- c("B1","B2","B3","B4","B5","B1","B2","B3","B4","B5","B1","B2","B3","B4","B5",
           "B1","B2","B3","B4","B5",
           "B1","B2","B3","B4","B5","B1","B2","B3","B4","B5","B1","B2","B3","B4","B5",
           "B1","B2","B3","B4","B5")
           
# Criando o dataframe e fazendo analise
dados_2 <- data.frame(Tratamento_2, Leite_2, Bloco_2)
anova4 <- aov(Leite_2 ~ Tratamento_2*Bloco_2, data = dados_2)
summary(anova4)




#Comparação entre as médias
TukeyHSD(anova)

duncan <- duncan.test(anova <- aov(Leite ~ Tratamento, data = dados), trt="Tratamento", group=T)
print(duncan$groups)
print(duncan$means)
ScheffeTest(anova)

# 1. Homogeneity of variances
plot(anova, 1)

leveneTest(Leite ~ Tratamento, data = dados)

bartlett.test(Leite ~ Tratamento, data = dados)

# Levene's test is an alternative to the Bartlett test. 

# H0 As variâncias são homogêneas
# H1 As variâncias não são homogêneas
Como p-valor calculado (p=0.3502)é maior que o nível de significância adotado
(α=0,05), não se rejeita H0. Logo, as variâncias são homogêneas.



# 2. Normality
shapiro.test(anova$res)
# H0: Os erros têm distribuição normal
# H1: Os erros não têm distribuição normal
# Como p-valor calculado (p=0.7564) é maior que o
# nível de significância 
# adotado (α=0,05 ), não se rejeita H0.
# Logo, os erros seguem distribuição normal.

plot(anova, 2)


###########  outros pacotes 
library("easyanova")
library("knitr")
library("ExpDes.pt")

with(dados,dic(dados$Tratamento, dados$Leite, mcomp="tukey"))
with(dados,dbc(dados$Tratamento, dados$Bloco,dados$Leite, mcomp="tukey"))
# opções para serem usadas no mcomp
# lsd, duncan , snk, bonferroni
### analisando a diferença em relação a média

mod1 = with(dados, aov(Leite ~ Tratamento))
mod2 = with(dados, aov(Leite ~ Tratamento+Bloco))

TukeyHSD(mod1, "Tratamento", ordered = TRUE)
TukeyHSD(mod2, "Tratamento", ordered = TRUE)

plot(TukeyHSD(mod1, "Tratamento"), col='blue', las=1)
plot(TukeyHSD(mod2, "Tratamento"), col='blue', las=1)

### DQL

Tratamento <- c("A", "A", "A", "A","B", "B", "B", "B","C", "C", "C", "C","D", "D", "D", "D")
Animal <- c(1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4)
Periodo <- c("I","III","IV","II","II","I","III","IV","IV","II","I","III","III","IV","II","I")
Leite <- c(23,21,22,20,25,22,23,20,26,27,20,23,27,28,25,24)

str(Leite)

dados <- data.frame(Tratamento, Animal, Leite, Periodo)
dados$Animal <- as.factor(dados$Animal)

mod <- aov(Leite ~ Animal + Periodo + Tratamento,dados) 
summary(mod)

# pacote ExpDes.pt

library("ExpDes.pt")

dql(
  trat=Tratamento,
  linha=Periodo,
  coluna=Animal,
  resp=Leite,
  mcomp="tukey",
  sigT= 0.05,
  sigF= 0.05,
  unfold=NULL
)


```

## Quadrado latino + analise de variancia

```{r}
# quadrado latino
library(agricolae)
library(gridExtra)
library(grid)


# Criando o croki de tratamentos
croqui=function(trat){
  r=length(trat)
  sort=design.lsd(trat,r,serie=0)
  sort$book[,4]=as.factor(matrix(sort$book[,4],r,,T))
  ncol=r
  gs <- lapply(sort$book[,4], function(ii)
    grobTree(rectGrob(gp=gpar(fill=ii, alpha=0.5)),textGrob(ii)))
  grid.arrange(grobs=gs, ncol=ncol)}

trat=c("T1","T2","T3","T4")

croqui(trat)

#Criando o dataframe

install.packages("kableExtra")
install.packages("hnp")
library(kableExtra)
library(hnp)

Conjunto de dados
RESP=c(93.0, 115.4, 116.9, 110.2, 110.4,110.6, 96.5, 108.9, 97.6, 112.0,102.1, 108.6, 77.9, 102.0, 111.7,115.4, 94.9, 114.0, 100.2, 118.5,117.6, 114.1, 118.7, 108.8, 80.2)
(TRAT=c("A","C","E","D","B","C","E","B","A","D","B","D","A","E","C","D","A","C","B","E","E","B","D","C","A"))
(linha=as.factor(rep(1:5,each=5)))
(coluna=as.factor(rep(1:5,5)))
dados = data.frame(TRAT, linha, coluna, RESP)
alfa=0.05

#Resumo dos dados
Media=mean(RESP)
Desvio=sd(RESP)
Variancia=var(RESP)
Maximo=max(RESP)
Minimo=min(RESP)
Mediana=median(RESP)
descritiva=cbind(Media,
                 Desvio, 
                 Variancia, 
                 Maximo, 
                 Minimo, 
                 Mediana)
kable(descritiva)

#Resumo dos dados por tratamento
Media=tapply(RESP,TRAT, mean)
Desvio=tapply(RESP,TRAT,sd)
Variancia=tapply(RESP,TRAT, var)
Maximo=tapply(RESP,TRAT,max)
Minimo=tapply(RESP,TRAT, min)
Mediana=tapply(RESP,TRAT,median)
descritiva=cbind(Media,
                 Desvio, 
                 Variancia, 
                 Maximo, 
                 Minimo, 
                 Mediana)
kable(descritiva)


# Fazendo o boxplot
car::Boxplot(RESP~TRAT,
             las=1,
             col="lightblue", xlab="",
             ylab=expression("Resposta"))

#Analise de variancia
points(Media,col="red", pch=8)


mod=aov(RESP~ TRAT+linha+coluna)
av=anova(mod)
names(av)=c("GL","SQ","QM","Teste F", "p-valor")
kable(av, align = "l", format="pandoc")

#suposicoes

#Normalidade
#####
# Ho: Os erros seguem distribuição normal
# H1: Os erros não seguem distribuição normal

#Shapiro wilk
(norm=shapiro.test(mod$res))

# Como p-valor calculado (p=0.438) é maior que o nível de significância adotado (α=0.05), não rejeita-se  
# Ho. Logo, os erros seguem distribuição normal..
hnp::hnp(mod, las=1, xlab="Quantis teóricos", pch=16)

#homogeneidade dos dados
(homog=with(dados, bartlett.test(mod$res ~ TRAT)))

# Ho: Os erros seguem distribuição normal
# H1: Os erros não seguem distribuição normal
# Como p-valor calculado (p=0.0996) é maior que o nível de significância adotado (p=0.05), não rejeita-se Ho. Logo, as variâncias são homogêneas.

#Independencia dos erros 
(ind=lmtest::dwtest(mod))
 
# H0: Os erros são independentes
# H1: Os erros nao são independentes

# distribuição dos erros 
plot(mod$res, las=1, pch=19, col='red')
abline(h=0)

# teste da anova
library(easyanova)

ea1(dados,design = 3)

#vendo o pacote
?ea1

# pacote laercio
install.packages("laercio")
require(laercio)
library(laercio)
LTukey(mod,"trat",conf.level=0.95)

#Pacte agricolae
require(agricolae)
TukeyHSD(mod, "TRAT", ordered = TRUE)

#Pacte ExpDes.pt
require(ExpDes.pt)
library(ExpDes.pt)
dql(TRAT,linha,coluna,RESP)
```


## Analise não paramétrica

```{r}
#### Aula – Análises Não Paramétricas em R ####
# Este script é um guia prático para os principais testes estatísticos não paramétricos em R.
# Todos os testes incluem os comparativos (post-hoc), mesmo quando o teste global não é significativo,
# para fins didáticos.

# Instalar e carregar pacotes necessários
# install.packages("rstatix")
# install.packages("ggplot2")
# install.packages("car")
# install.packages("ggpubr")
# install.packages("dplyr")
# install.packages("tidyr") # Para organizar dados para testes pareados
library(rstatix)
library(ggplot2)
library(car)
library(ggpubr)
library(dplyr)
library(tidyr)

### 1. Associação entre variáveis categóricas
# Qui-quadrado + Fisher
raca <- c(rep("Angus", 25), rep("Meio sangue", 35), rep("Nelore", 25))
prenhes <- c(rep(1, 10), rep(0, 15),
             rep(1, 4), rep(0, 31),
             rep(1, 8), rep(0, 17))
dataframe <- data.frame(Raca = as.factor(raca), Prenhes = as.factor(prenhes))

tabela <- table(dataframe$Raca, dataframe$Prenhes)
print(tabela)

barplot(tabela, beside = TRUE, legend.text = rownames(tabela),
        main = "Contagem de Prenhez por Raça", xlab = "Status de Prenhez (0=Não, 1=Sim)",
        col = c("lightblue", "orange", "lightgreen"))

resultados_chi2 <- chisq.test(tabela)
print(resultados_chi2)

resultados_fisher <- fisher.test(tabela)
print(resultados_fisher)


### 2. Dois grupos independentes
# Mann–Whitney (Wilcoxon rank-sum)
grupo <- c(rep("Fertilizante A", 8), rep("Fertilizante B", 8))
altura <- c(21, 22, 20, 19, 23, 22, 21, 20,
            25, 24, 23, 22, 26, 27, 25, 24)
dados_mannwhitney <- data.frame(grupo, altura)

ggplot(dados_mannwhitney, aes(x = grupo, y = altura, fill = grupo)) +
  geom_boxplot() +
  labs(title = "Altura das Plantas por Tipo de Fertilizante", y = "Altura (cm)") +
  theme_minimal() +
  stat_compare_means(method = "wilcox.test", label.y = 28)

resultados_mannwhitney <- wilcox.test(altura ~ grupo, data = dados_mannwhitney)
print(resultados_mannwhitney)

# Comparação par a par (mesmo com 2 grupos, mostra p-valor ajustado)
pairwise_result <- pairwise.wilcox.test(dados_mannwhitney$altura, dados_mannwhitney$grupo,
                                        p.adjust.method = "bonferroni")
print(pairwise_result)


### 3. Dois grupos dependentes (pareados)
# Wilcoxon signed-rank
antes <- c(80, 75, 90, 85, 95, 78, 82, 88)
depois <- c(78, 73, 88, 83, 92, 75, 80, 85)
dados_wilcoxon_pareado <- data.frame(antes, depois)

resultados_wilcoxon_pareado <- wilcox.test(antes, depois, paired = TRUE)
print(resultados_wilcoxon_pareado)

# Comparação par a par (forçado para fins didáticos)
pairwise_pareado <- pairwise_wilcox_test(
  data = tidyr::pivot_longer(dados_wilcoxon_pareado, cols = c(antes, depois),
                             names_to = "tempo", values_to = "peso"),
  peso ~ tempo, paired = TRUE, p.adjust.method = "bonferroni"
)
print(pairwise_pareado)


### 4. Mais de dois grupos independentes
# Kruskal–Wallis + Dunn
Dieta <- c(rep("A", 5), rep("B", 5), rep("C", 5), rep("D", 5))
Producao <- c(23, 21, 22, 20, 21,
              25, 22, 23, 20, 19,
              26, 27, 20, 24, 22,
              27, 28, 25, 23, 24)
dados_kruskal <- data.frame(Dieta, Producao)

ggplot(dados_kruskal, aes(x = Dieta, y = Producao, fill = Dieta)) +
  geom_boxplot() +
  labs(title = "Produção de Leite por Dieta") +
  theme_minimal()

kruskal_result <- kruskal.test(Producao ~ Dieta, data = dados_kruskal)
print(kruskal_result)

posthoc_dunn <- dunn_test(Producao ~ Dieta, data = dados_kruskal, p.adjust.method = "bonferroni")
print(posthoc_dunn %>% select(group1, group2, p.adj, p.adj.signif))


### 5. Mais de dois grupos dependentes
# Friedman + Wilcoxon pareado post-hoc
vaca <- factor(rep(1:6, each = 3))
dieta <- factor(rep(c("A", "B", "C"), times = 6))
producao <- c(20, 22, 24, 18, 20, 19,
              23, 25, 26, 19, 18, 20,
              21, 23, 24, 20, 22, 23)
dados_friedman <- data.frame(vaca, dieta, producao)

friedman_result <- friedman.test(producao ~ dieta | vaca, data = dados_friedman)
print(friedman_result)

posthoc_friedman <- pairwise_wilcox_test(producao ~ dieta, data = dados_friedman,
                                         paired = TRUE, p.adjust.method = "bonferroni")
print(posthoc_friedman %>% select(group1, group2, p.adj, p.adj.signif))


### 6. Testes de Pressupostos (Normalidade e Homogeneidade)
# Este trecho de código assume que você já rodou um modelo linear
# com a função `aov` ou `lm`. Por exemplo: `modelo <- aov(peso ~ grupo, data = dados)`

# Instale o pacote `car` se ainda não tiver
# install.packages("car")

# Exemplo de dados para o teste de Levene e Bartlett
grupo_levene <- rep(c("A", "B", "C"), each = 5)
peso <- c(21, 22, 20, 19, 23,
          25, 24, 23, 22, 26,
          27, 28, 26, 25, 29)
dados_levene <- data.frame(grupo = grupo_levene, peso = peso)

# Exemplo de modelo para extrair os resíduos
modelo_exemplo <- aov(peso ~ grupo, data = dados_levene)

# Testes de Normalidade dos Resíduos
# Teste de Shapiro-Wilk (para amostras até 50)
shapiro_result <- shapiro.test(resid(modelo_exemplo))
print(shapiro_result)

# Teste de Kolmogorov-Smirnov (para amostras > 50)
ks_result <- ks.test(resid(modelo_exemplo), "pnorm")
print(ks_result)

# Testes de Homogeneidade (Homocedasticidade)
# Teste de Levene (robusto a não-normalidade)
levene_result <- leveneTest(peso ~ grupo, data = dados_levene)
print(levene_result)

# Teste de Bartlett (para dados com distribuição normal)
bartlett_result <- bartlett.test(peso ~ grupo, data = dados_levene)
print(bartlett_result)
```



# Analise de correlação e correspondência

[script](https://drive.google.com/drive/folders/1q6dx9VHyB-resROL8dXAymwLQG4Hvcyh?usp=sharing)
Aula

```{r}
# ==============================================================================
# ROTEIRO DE AULA: ANÁLISE DE CORRELAÇÃO, COM GRÁFICOS E HEATMAPS
# ==============================================================================

# Pacotes necessários
pacotes <- c("car","ggpubr","corrplot","PerformanceAnalytics","GGally",
             "psych","qgraph","gplots","patchwork","pheatmap")

# Função para checar disponibilidade dos pacotes
check_packages_availability <- function(packages) {
  availability <- sapply(packages, requireNamespace, quietly = TRUE)
  data.frame(Pacote = names(availability), Disponivel = unname(availability))
}
check_packages_availability(pacotes)

# Carregar pacotes
lapply(pacotes, library, character.only = TRUE)

# =========================
# 1. Análise com a base mtcars
# =========================
my_data <- mtcars
cat("Visualizando as primeiras linhas:\n")
head(my_data, 6)

# 1.1 Gráfico de dispersão mpg x wt
g1 <- ggscatter(my_data, x = "mpg", y = "wt", 
                add = "reg.line", conf.int = TRUE, 
                cor.coef = TRUE, cor.method = "pearson",
                xlab = "Consumo (milhas/galão)", 
                ylab = "Peso (1000 lbs)",
                title = "mpg vs wt - Correlação Pearson")

# 1.2 Testes de normalidade Shapiro-Wilk e Q-Q plots
shapiro_mpg <- shapiro.test(my_data$mpg)
shapiro_wt <- shapiro.test(my_data$wt)
cat("Shapiro mpg:", shapiro_mpg$p.value, "\n")
cat("Shapiro wt:", shapiro_wt$p.value, "\n")

gqq_mpg <- ggqqplot(my_data$mpg, ylab="MPG")
gqq_wt  <- ggqqplot(my_data$wt, ylab="WT")
gqq_mpg + gqq_wt + plot_annotation(title="Q-Q Plots para Normalidade")

# 1.3 Correlação entre mpg e wt
res_pearson  <- cor.test(my_data$wt, my_data$mpg, method="pearson")
res_spearman <- cor.test(my_data$wt, my_data$mpg, method="spearman")
cat("Pearson:", res_pearson$estimate, "p =", res_pearson$p.value, "\n")
cat("Spearman:", res_spearman$estimate, "p =", res_spearman$p.value, "\n")

# 1.4 Matriz de correlação e visualizações
vars <- c("mpg","cyl","disp","hp","wt","qsec")
sub_data <- my_data[, vars]
M_corr <- cor(sub_data)

# Corrplots com valores
corrplot(M_corr, method="circle", type="upper", addCoef.col="black", number.cex=0.8,
         tl.col="black", tl.srt=45, title="Matriz de Correlação (Pearson)")

# Correlação 2 a 2 com histogramas
chart.Correlation(sub_data, histogram=TRUE, pch=19)

# Heatmaps com pheatmap
col_palette <- colorRampPalette(c("blue","white","red"))(50)
pheatmap(M_corr, color=col_palette, display_numbers=TRUE, number_format="%.2f",
         cluster_rows=TRUE, cluster_cols=TRUE, main="Heatmap Correlação (Pearson)")
pheatmap(scale(sub_data), color=col_palette, display_numbers=TRUE, number_format="%.2f",
         cluster_rows=TRUE, cluster_cols=TRUE, main="Heatmap Dados Escalonados")

# =========================
# 2. Análise com dados de agronomia
# =========================
DPF=c(46,46,46,46,46,46,43,43,43,46)
APF=c(58.33,55,50,41,35.67,43.33,35.67,36,35.33,46.67)
DPM=c(105,105,102,110,110,112,110,110,105,112)
APM=c(100,90.33,97,91.33,97.67,77.33,90,93,91.33,98)
IPV=c(15,20,17,10,22.67,14.33,23,19.33,15.33,14.33)
ACA=c(2,1.9,2.2,1.5,1.2,1,2,1.5,1.2,3)
PRO=c(2444.44,2870.37,2314.81,2629.63,2444.44,2592.59,2962.96,3037.04,3037.04,2592.59)
MCG=c(10.78,10.96,10.07,10.77,11.17,11.24,12.57,13.35,13.77,14.23)

dados <- data.frame(DPF,APF,DPM,APM,IPV,ACA,PRO,MCG)
head(dados)

# Matrizes Pearson e Spearman
M_pearson  <- cor(dados, method="pearson")
M_spearman <- cor(dados, method="spearman")

# Função para corrplot com valores
plot_corr <- function(M, method_name){
  corrplot(M, method="circle", type="upper", addCoef.col="black", number.cex=0.8,
           tl.col="black", tl.srt=45, title=paste("Correlação", method_name))
}

plot_corr(M_pearson, "Pearson")
plot_corr(M_spearman, "Spearman")

# Pairs.panels
pairs.panels(dados, method="pearson", hist.col="blue", density=TRUE, ellipses=TRUE, stars=TRUE)
pairs.panels(dados, method="spearman", hist.col="red", density=TRUE, ellipses=TRUE, stars=TRUE)

# Redes de correlação com qgraph
qgraph(M_pearson, shape="circle", posCol="darkgreen", negCol="darkred", layout="groups", vsize=10)
qgraph(M_spearman, shape="circle", posCol="darkgreen", negCol="darkred", layout="groups", vsize=10)

```




# Analise de Cluster
[Aula]()
[Rscript](https://drive.google.com/file/d/1peP7Uz85eiZFWJdHxdiCWzvog7ATHzJ6/view?usp=sharing)

```{r}
########################################
#
#   CHAMANDO BIBLIOTECAS IMPORTANTES
#
########################################

library(tidyverse) #pacote para manipulacao de dados
library(cluster) #algoritmo de cluster
library(dendextend) #compara dendogramas
library(factoextra) #algoritmo de cluster e visualizacao
library(fpc) #algoritmo de cluster e visualizacao
library(gridExtra) #para a funcao grid arrange
library(readxl)


########################################
#
#         CLUSTER HIERARQUICO - juntos
#
########################################

#LEITURA DOS DADOS
alunos_pap <- read.table("dados/alunos_pap.csv", sep = ";", header = T, dec = ",")
View(alunos_pap)
rownames(alunos_pap) <- alunos_pap[,1]
alunos_pap <- alunos_pap[,-1]

#CALCULANDO MATRIZ DE DISTANCIAS
d <- dist(alunos_pap, method = "euclidean")
d

#DEFININDO O CLUSTER A PARTIR DO METODO ESCOLHIDO
#metodos disponiveis "average", "single", "complete" e "ward.D"
hc1 <- hclust(d, method = "single" )
hc2 <- hclust(d, method = "complete" )
hc3 <- hclust(d, method = "average" )
hc4 <- hclust(d, method = "ward.D" )

#DESENHANDO O DENDOGRAMA
plot(hc1, cex = 0.6, hang = -1)
plot(hc2, cex = 0.6, hang = -1)
plot(hc3, cex = 0.6, hang = -1)
plot(hc4, cex = 0.6, hang = -1)

#BRINCANDO COM O DENDOGRAMA PARA 2 GRUPOS
rect.hclust(hc4, k = 2)

#COMPARANDO DENDOGRAMAS
#comparando o metodo average com ward
dend3 <- as.dendrogram(hc3)
dend4 <- as.dendrogram(hc4)
dend_list <- dendlist(dend3, dend4) 
#EMARANHADO, quanto menor, mais iguais os dendogramas sao
tanglegram(dend3, dend4, main = paste("Emaranhado =", round(entanglement(dend_list),2)))
#agora comparando o metodo single com complete
dend1 <- as.dendrogram(hc1)
dend2 <- as.dendrogram(hc2)
dend_list2 <- dendlist(dend1, dend2) 
#EMARANHADO, quanto menor, mais iguais os dendogramas sao
tanglegram(dend1, dend2, main = paste("Emaranhado =", round(entanglement(dend_list2),2)))

#criando 2 grupos de alunos
grupo_alunos2 <- cutree(hc4, k = 2)
table(grupo_alunos2)

#transformando em data frame a saida do cluster
alunos_grupos <- data.frame(grupo_alunos2)

#juntando com a base original
Base_alunos_fim <- cbind(alunos_pap, alunos_grupos)

# entendendo os clusters
#FAZENDO ANALISE DESCRITIVA
#MEDIAS das variaveis por grupo
mediagrupo_alunos <- Base_alunos_fim %>% 
  group_by(grupo_alunos2) %>% 
  summarise(n = n(),
            Portugues = mean(Portugues), 
            Matematica = mean(Matematica))
mediagrupo_alunos


########################################
#
#     CLUSTER HIERARQUICO - MCDonald
#
########################################

#Carregar base de dados: 
mcdonalds <- read.table("dados/MCDONALDS.csv", sep = ";", dec = ",", header = T)
#transformar o nome dos lanches em linhas
rownames(mcdonalds) <- mcdonalds[,1]
mcdonalds <- mcdonalds[,-1]

#Padronizar variaveis
mcdonalds.padronizado <- scale(mcdonalds)

#calcular as distancias da matriz utilizando a distancia euclidiana
distancia <- dist(mcdonalds.padronizado, method = "euclidean")

#Calcular o Cluster: metodos disponiveis "average", "single", "complete" e "ward.D"
cluster.hierarquico <- hclust(distancia, method = "single" )

# Dendrograma
plot(cluster.hierarquico, cex = 0.6, hang = -1)

#Criar o grafico e destacar os grupos
rect.hclust(cluster.hierarquico, k = 2)

#VERIFICANDO ELBOW 
fviz_nbclust(mcdonalds.padronizado, FUN = hcut, method = "wss")


#criando 4 grupos de lanches
grupo_lanches4 <- cutree(cluster.hierarquico, k = 4)
table(grupo_lanches4)

#transformando em data frame a saida do cluster
Lanches_grupos <- data.frame(grupo_lanches4)

#juntando com a base original
Base_lanches_fim <- cbind(mcdonalds, Lanches_grupos)

#FAZENDO ANALISE DESCRITIVA
#MEDIAS das variaveis por grupo
mediagrupo <- Base_lanches_fim %>% 
  group_by(grupo_lanches4) %>% 
  summarise(n = n(),
            Valor.Energetico = mean(Valor.Energetico), 
            Carboidratos = mean(Carboidratos), 
            Proteinas = mean(Proteinas),
            Gorduras.Totais = mean(Gorduras.Totais), 
            Gorduras.Saturadas = mean(Gorduras.Saturadas), 
            Gorduras.Trans = mean(Gorduras.Trans),
            Colesterol = mean(Colesterol), 
            Fibra.Alimentar = mean(Fibra.Alimentar), 
            Sodio = mean(Sodio),
            Calcio = mean(Calcio), 
            Ferro = mean(Ferro) )
mediagrupo


#    CLUSTER NAO HIERARQUICO - Mcdonald


#AGRUPANDO LANCHES PELO METODO NAO HIERARQUICO

#Rodar o modelo
mcdonalds.k2 <- kmeans(mcdonalds.padronizado, centers = 2)

#Visualizar os clusters
fviz_cluster(mcdonalds.k2, data = mcdonalds.padronizado, main = "Cluster K2")

#Criar clusters
mcdonalds.k3 <- kmeans(mcdonalds.padronizado, centers = 3)
mcdonalds.k4 <- kmeans(mcdonalds.padronizado, centers = 4)
mcdonalds.k5 <- kmeans(mcdonalds.padronizado, centers = 5)

#Criar graficos
G1 <- fviz_cluster(mcdonalds.k2, geom = "point", data = mcdonalds.padronizado) + ggtitle("k = 2")
G2 <- fviz_cluster(mcdonalds.k3, geom = "point",  data = mcdonalds.padronizado) + ggtitle("k = 3")
G3 <- fviz_cluster(mcdonalds.k4, geom = "point",  data = mcdonalds.padronizado) + ggtitle("k = 4")
G4 <- fviz_cluster(mcdonalds.k5, geom = "point",  data = mcdonalds.padronizado) + ggtitle("k = 5")

#Imprimir graficos na mesma tela
grid.arrange(G1, G2, G3, G4, nrow = 2)

#VERIFICANDO ELBOW 
fviz_nbclust(mcdonalds.padronizado, kmeans, method = "wss")


########################################
#
#    CLUSTER NAO HIERARQUICO - Municipios
#
########################################

#carregar base municipio
municipios <- read.table("dados/municipios.csv", sep = ";", header = T, dec = ",")
rownames(municipios) <- municipios[,1]
municipios <- municipios[,-1]


#padronizar dados
municipios.padronizado <- scale(municipios)


#Agora vamos rodar de 3 a 6 centros  e visualizar qual a melhor divisao
municipios.k3 <- kmeans(municipios.padronizado, centers = 3)
municipios.k4 <- kmeans(municipios.padronizado, centers = 4)
municipios.k5 <- kmeans(municipios.padronizado, centers = 5)
municipios.k6 <- kmeans(municipios.padronizado, centers = 6)

#Graficos
G1 <- fviz_cluster(municipios.k3, geom = "point", data = municipios.padronizado) + ggtitle("k = 3")
G2 <- fviz_cluster(municipios.k4, geom = "point",  data = municipios.padronizado) + ggtitle("k = 4")
G3 <- fviz_cluster(municipios.k5, geom = "point",  data = municipios.padronizado) + ggtitle("k = 5")
G4 <- fviz_cluster(municipios.k6, geom = "point",  data = municipios.padronizado) + ggtitle("k = 6")

#Criar uma matriz com 4 graficos
grid.arrange(G1, G2, G3, G4, nrow = 2)

#VERIFICANDO ELBOW 
fviz_nbclust(municipios.padronizado, FUN = hcut, method = "wss")

#juntando dados
municipios2 <- read.table("dados/municipios.csv", sep = ";", header = T, dec = ",")
municipiosfit <- data.frame(municipios.k6$cluster)

#Agrupar cluster e base
MunicipioFinal <-  cbind(municipios2, municipiosfit)

```

# Analise de componentes principais PCA
[Teorica e Prática](https://drive.google.com/drive/folders/18cxLyCvyyIWKyEIwKWRI3AlWgQUkywV-?usp=sharing)

# Regressão linear - Simples e Múltipla
Regressõa simples

```{r}
###############################################################

# ============================================================
# 1. Introdução teórica (explicação para os alunos)
# ============================================================

# A regressão linear simples é usada quando queremos analisar a
# relação entre duas variáveis quantitativas:
#   - Uma variável independente (X): explicativa, preditora
#   - Uma variável dependente (Y): resposta, desfecho

# Exemplo veterinário:
#   - X: Idade do animal (em meses)
#   - Y: Peso corporal (em kg)
#   Pergunta: "O peso aumenta linearmente com a idade?"

# ============================================================
# 2. Criando o conjunto de dados
# ============================================================


# Limpando o ambiente
rm(list = ls())
cat("\014")

# Desativar notação científica
options(scipen = 999)


# Dados simulados representando 15 bezerros
# (idade em meses e peso em kg)
dados <- data.frame(
  idade_meses = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15),
  peso_kg = c(38, 45, 52, 60, 65, 72, 80, 84, 90, 96, 102, 106, 112, 117, 122)
)

# Visualizar os primeiros dados
print(dados)

# ============================================================
# 3. Análise exploratória
# ============================================================

# Resumo descritivo
summary(dados)

# Gráfico exploratório (dispersão)
plot(dados$idade_meses, dados$peso_kg,
     main = "Relação entre Idade e Peso de Bezerros",
     xlab = "Idade (meses)",
     ylab = "Peso corporal (kg)",
     pch = 19, col = "blue")

# Comentário para o aluno:
# Observem que os pontos parecem seguir uma tendência crescente.
# Isso sugere uma relação linear positiva: quanto maior a idade, maior o peso.

# ============================================================
# 4. Ajustando o modelo de regressão linear simples
# ============================================================

modelo <- lm(peso_kg ~ idade_meses, data = dados)

# Ver o resumo completo do modelo
summary(modelo)

# Interpretação dos principais resultados:
# - Coeficiente angular (β1): indica o quanto o peso aumenta, em média,
#   a cada 1 mês de idade.
# - Intercepto (β0): peso estimado quando a idade = 0 (geralmente interpretado
#   apenas matematicamente).
# - R²: porcentagem da variação do peso explicada pela idade.
# - p-valor: testa se a relação entre idade e peso é estatisticamente significativa.

# ============================================================
# 5. Visualizando a reta de regressão
# ============================================================

# Adicionando a reta de regressão ao gráfico
abline(modelo, col = "red", lwd = 2)

# Legenda
legend("topleft",
       legend = c("Dados observados", "Reta ajustada"),
       col = c("blue", "red"),
       pch = c(19, NA),
       lty = c(NA, 1),
       bty = "n")

# ============================================================
# 6. Interpretando o modelo numericamente
# ============================================================

coef(modelo)  # mostra os coeficientes (intercepto e inclinação)

# Podemos escrever a equação do modelo:
# Peso = β0 + β1 * Idade
# Substituindo pelos valores obtidos:
# Exemplo: Peso = 35 + 6.2 * Idade

# Isso significa que a cada mês, o peso médio aumenta 6.2 kg.

# ============================================================
# 7. Fazendo previsões
# ============================================================

# Vamos prever o peso de bezerros com 8, 12 e 16 meses:
novos_dados <- data.frame(idade_meses = c(8, 12, 16))
predicoes <- predict(modelo, novos_dados, interval = "confidence")
predicoes

# Intervalos de confiança ajudam a entender a incerteza da previsão.

# ============================================================
# 8. Visualizando as previsões
# ============================================================

# Vamos adicionar as previsões ao gráfico
plot(dados$idade_meses, dados$peso_kg,
     main = "Regressão Linear Simples - Bezerros",
     xlab = "Idade (meses)",
     ylab = "Peso corporal (kg)",
     pch = 19, col = "blue")

# Reta ajustada
abline(modelo, col = "red", lwd = 2)

# Intervalos de confiança
ic <- predict(modelo, interval = "confidence")
lines(dados$idade_meses, ic[, "lwr"], col = "gray", lty = 2)
lines(dados$idade_meses, ic[, "upr"], col = "gray", lty = 2)

# Comentário:
# As linhas cinzas representam os limites superior e inferior do intervalo de confiança
# para a média estimada — mostram a incerteza na previsão do peso médio.

# ============================================================
# 9. Diagnóstico do modelo
# ============================================================

# Testar os pressupostos da regressão:
# - Resíduos devem ter distribuição aproximadamente normal
# - Homocedasticidade (variância constante)
# - Relação linear entre as variáveis

par(mfrow = c(2, 2))
plot(modelo)
par(mfrow = c(1, 1))

# Comentário:
# Se os resíduos estiverem aleatoriamente espalhados (sem padrão),
# e o gráfico "Normal Q-Q" mostrar uma linha aproximadamente reta,
# o modelo linear é adequado.

#Residuals vs fitted - Curva em U
# Objetivo: verificar linearidade e homocedasticidade (variância constante dos erros).
#Os pontos devem estar espalhados aleatoriamente em torno da linha zero (sem forma definida).
#Isso indica que o modelo linear explica bem os dados.
# a idade pode não ser completamente linear  
# crescimento pode desacelerar em idade maiores

#Normal Q-
#Objetivo: verificar normalidade dos resíduos — outro pressuposto da regressão.
#Espeado Os pontos devem cair aproximadamente sobre a linha reta tracejada.
#Os pontos estão quase sobre a linha, mas com pequenas curvaturas nas extremidades
# indicando ligeiro desvio da normalidade, mas nada grave.

#Scale location
# Objetivo: verificar se os resíduos têm variância constante (homocedasticidade).
# Esperado : Pontos aleatórios espalhados sem padrão e sem aumento da dispersão.
#A linha vermelha está ligeiramente curvada, sugerindo heterocedasticidade leve 
# os resíduos aumentam ou diminuem um pouco conforme o valor previsto.

#Residuals vs Leverage
# Objetivo: identificar pontos influentes — observações que exercem forte impacto sobre a reta ajustada.
#O que esperamos ver:
#  A maioria dos pontos deve estar próxima do centro, sem ultrapassar as linhas tracejadas de Cook’s distance.
#O ponto 15 (idade mais alta) parece ter alta alavancagem (distante dos demais), mas não ultrapassa a linha de Cook
# portanto, não é um ponto influente grave, apenas merece atenção.


##### Distância de Cook
# A Cook’s distance (ou distância de Cook) 
# é uma medida que mostra quanto cada ponto individual influencia o ajuste do modelo de regressão.

# ============================================================
# 10. Avaliando a qualidade do ajuste
# ============================================================

r2 <- summary(modelo)$r.squared
cat("Coeficiente de determinação (R²):", round(r2, 3), "\n")

# Interpretação:
# Um R² de 0.95 significa que 95% da variação no peso
# é explicada pela idade do bezerro.

############## Exercicio 


# ==========================================================
# Exercício 2: Regressão Linear Simples
# Tema: Produção de leite em função dos dias de lactação
# ==========================================================

# Criação do conjunto de dados
dados_leite <- data.frame(
  vaca = 1:15,
  dias_lactacao = c(10, 30, 50, 70, 90, 110, 130, 150, 170, 190, 210, 230, 250, 270, 290),
  leite_L = c(18, 22, 26, 29, 32, 33, 34, 33, 31, 28, 25, 22, 18, 14, 10) # curva típica de lactação
)

# Visualizar os dados
print(dados_leite)

# 1. Gráfico exploratório
plot(dados_leite$dias_lactacao, dados_leite$leite_L,
     main = "Produção de Leite vs Dias de Lactação",
     xlab = "Dias de Lactação",
     ylab = "Leite (L/dia)",
     pch = 19, col = "darkgreen")

# 2. Ajuste do modelo de regressão linear
modelo_leite <- lm(leite_L ~ dias_lactacao, data = dados_leite)

# 3. Resumo do modelo
summary(modelo_leite)

# 4. Adicionar a reta de regressão
abline(modelo_leite, col = "red", lwd = 2)

# 5. Diagnóstico do modelo
par(mfrow = c(2, 2))
plot(modelo_leite)

# 6. Cook’s Distance
cooks_leite <- cooks.distance(modelo_leite)
print(cooks_leite)

# Critério prático (regra geral)
# 
# 🔹 Cook’s Distance < 0.5 → ponto não influente, seguro.
# 
# 🟡 Entre 0.5 e 1.0 → ponto moderadamente influente, precisa de atenção.
# 
# 🔴 > 1.0 → ponto muito influente, pode estar distorcendo o modelo.

plot(cooks_leite, type = "h", main = "Cook’s Distance - Produção de Leite", ylab = "Influência", col = "darkred")
abline(h = 0.5, col = "blue", lty = 2)

# 7. Intervalos de confiança
ic_leite <- predict(modelo_leite, interval = "confidence")
head(ic_leite)

r2_leite <- summary(modelo_leite)$r.squared
cat("Coeficiente de determinação (R²):", round(r2_leite, 3), "\n")


# 8. Perguntas para análise
# - A relação entre dias de lactação e produção de leite parece linear?
# - O modelo explica bem a variação nos dados (veja o R²)?
# - Há curvatura nos resíduos, indicando relação não-linear?
# - Algum ponto é influente segundo a Cook’s Distance?
# - Como você interpretaria o coeficiente angular neste contexto?

####
# ============================================================
# 11. Teste de Breusch-Pagan (Heterocedasticidade)
# ============================================================

# O teste de Breusch-Pagan avalia se a variância dos resíduos
# é constante (homocedasticidade) ou não (heterocedasticidade).

# Pressuposto da regressão linear:
# A variância dos erros deve ser constante ao longo dos valores ajustados.
# Quando isso não ocorre, dizemos que há heterocedasticidade — 
# o que pode tornar os testes t e F inválidos (erro padrão incorreto).

# Para realizar o teste, precisamos carregar o pacote 'lmtest'
if(!require(lmtest)) install.packages("lmtest")
library(lmtest)

# Aplicando o teste ao modelo
bptest(modelo)

# Interpretação:
# - H0 (hipótese nula): os resíduos têm variância constante (homocedasticidade)
# - H1 (hipótese alternativa): há heterocedasticidade (variância não constante)
#
# Se o p-valor > 0.05 → NÃO rejeitamos H0 → variância é constante → OK 👍
# Se o p-valor < 0.05 → rejeitamos H0 → há heterocedasticidade → ⚠️ problema

# Exemplo de saída típica:
# studentized Breusch-Pagan test
# data:  modelo
# BP = 1.56, df = 1, p-value = 0.21
#
# Interpretação:
# Como p = 0.21 > 0.05, não há evidências de heterocedasticidade.
# Logo, o pressuposto de homocedasticidade é atendido.

# ============================================================
# 12. Aplicando o teste no segundo exemplo (produção de leite)
# ============================================================

bptest(modelo_leite)

# Exemplo de interpretação:
# Se o p-valor for < 0.05, significa que a variância dos resíduos muda
# ao longo dos dias de lactação — ou seja, o modelo linear pode não
# ser adequado (pode haver curvatura).
#
# 
# “Como a produção de leite segue uma curva típica (sobe e depois desce),
# é esperado que o modelo linear simples não capture bem a variação — 
# e o teste de Breusch-Pagan pode indicar heterocedasticidade.”

# ============================================================
# Conclusão
# ============================================================

# 🔍 Resumo do uso do Breusch-Pagan:
# - É um teste para confirmar o que observamos visualmente no gráfico
#   "Scale-Location" dos diagnósticos.
# - Serve para reforçar a análise estatística do pressuposto de variância constante.
# - Não substitui o gráfico, mas complementa a avaliação.

# 💡 Dica prática:
# Se houver heterocedasticidade (p < 0.05), podemos tentar:
# 1. Transformar a variável resposta (ex: log(Y))
# 2. Usar modelos robustos (ex: regressão robusta ou ponderada)
# 3. Ajustar um modelo não-linear (ex: quadrático ou polinomial)

```

Regressão linear múltipla

```{r}
###############################################################
# ============================================================
# 1. Preparação do ambiente
# ============================================================
rm(list = ls())
cat("\014")
options(scipen = 999)

# ============================================================
# 2. Criando o conjunto de dados (simulado)
# ============================================================
dados <- data.frame(
  idade_meses = 1:15,
  comprimento_cm = c(40,42,45,48,50,52,55,58,60,63,65,68,70,73,75),
  sexo = c(0,1,0,1,0,1,0,1,0,1,0,1,0,1,0),
  peso_kg = c(3.2,4.0,3.8,4.5,4.2,5.0,4.8,5.5,5.2,6.0,5.8,6.5,6.2,7.0,6.8)
)

# Visualizando dados
head(dados)
summary(dados)

# ============================================================
# 3. Gráficos exploratórios avançados
# ============================================================
library(ggplot2)

# Peso x Idade, colorindo por sexo
ggplot(dados, aes(x=idade_meses, y=peso_kg, color=factor(sexo))) +
  geom_point(size=3) +
  geom_smooth(method="lm", se=TRUE) +
  labs(x="Idade (meses)", y="Peso (kg)", color="Sexo") +
  ggtitle("Peso vs Idade por Sexo") +
  theme_minimal()

# Peso x Comprimento, colorindo por sexo
ggplot(dados, aes(x=comprimento_cm, y=peso_kg, color=factor(sexo))) +
  geom_point(size=3) +
  geom_smooth(method="lm", se=TRUE) +
  labs(x="Comprimento (cm)", y="Peso (kg)", color="Sexo") +
  ggtitle("Peso vs Comprimento por Sexo") +
  theme_minimal()

# ============================================================
# 4. Ajustando o modelo de regressão linear múltipla
# ============================================================
modelo <- lm(peso_kg ~ idade_meses + comprimento_cm + sexo, data=dados)
summary(modelo)

# ============================================================
# 5. Interpretando coeficientes
# ============================================================
coef(modelo)           # coeficientes do modelo
confint(modelo)        # intervalos de confiança dos coeficientes

#Interpretação do 2,5% e 97,5%:
# 2,5%: limite inferior do intervalo — o menor valor plausível para o coeficiente com 95% de confiança.
# 97,5%: limite superior do intervalo — o maior valor plausível para o coeficiente com 95% de confiança.

# #Exemplo com idade:
# 0.0625   0.4885
# O coeficiente verdadeiro do efeito da idade no peso provavelmente está entre 0,0625 e 0,4885 kg por mês, com 95% de confiança.
# Como o intervalo não inclui zero, concluímos que o efeito é estatisticamente significativo.
# Se o intervalo incluísse zero, não teríamos evidência suficiente de que a variável realmente afeta o peso.
  

# ============================================================
# 6. Diagnóstico do modelo
# ============================================================
par(mfrow=c(2,2))
plot(modelo)
par(mfrow=c(1,1))

# Cook's Distance
cooks <- cooks.distance(modelo)
plot(cooks, type="h", main="Cook's Distance", ylab="Influência", col="red")
abline(h=0.5, col="blue", lty=2)

# Perguntas:
# 4) Algum ponto é altamente influente?
# 5) O modelo atende às suposições de regressão linear?

# | Gráfico                                     | Objetivo                                  | Observação                                                                     | Interpretação                                        |
# | ------------------------------------------- | ----------------------------------------- | ------------------------------------------------------------------------------ | ---------------------------------------------------- |
# | **Residuals vs Fitted**                     | Verificar linearidade e homocedasticidade | Pontos dispersos aleatoriamente em torno da linha horizontal, sem padrão claro | Linearidade aceitável, sem forte heterocedasticidade |
# | **Q-Q Residuals**                           | Verificar normalidade dos resíduos        | A maioria dos pontos segue a linha teórica; pontos 3 e 7 levemente afastados   | Resíduos aproximadamente normais                     |
# | **Scale-Location**                          | Verificar homocedasticidade               | Não há padrão claro de aumento/diminuição da variância                         | Homocedasticidade aproximadamente respeitada         |
# | **Residuals vs Leverage (Cook’s Distance)** | Identificar pontos influentes             | Pontos 3 e 13 com maior leverage, mas Cook < 1                                 | Nenhum ponto crítico; atenção aos pontos 3, 7 e 13   |


# ============================================================
# 7. Previsões com intervalo de confiança
# ============================================================
novos <- data.frame(
  idade_meses=c(6,10,14),
  comprimento_cm=c(51,63,73),
  sexo=c(0,1,1)
)
predicoes <- predict(modelo, novos, interval="confidence")
predicoes

# Perguntas:
# 6) Qual a faixa de confiança das previsões?
# 7) Como interpretar o intervalo de confiança?

# ============================================================
# 8. Exercícios 
# ============================================================
# a) Criar um gráfico 3D ou scatter plot matrix com todas as variáveis
library(GGally)
ggpairs(dados, aes(color=factor(sexo)))

# b) Ajustar um modelo sem a variável sexo e comparar R²
modelo2 <- lm(peso_kg ~ idade_meses + comprimento_cm, data=dados)
summary(modelo2)

# Perguntas:
# 8) O R² mudou? O que isso significa sobre o efeito do sexo?
# 9) As variáveis idade e comprimento ainda são significativas?

###
# Ao retirar "sexo", idade e comprimento deixaram de ser significativos, pois os valores-p ficaram acima de 0,05. 
# Isso sugere que, sem sexo na equação, nenhuma das duas variáveis isoladamente explica estatisticamente o peso dos animais.
# Quando "sexo" está presente, ele absorve boa parte da explicação; 
# a força do efeito de idade aparece, enquanto comprimento permanece pouco relevante

####### adicional

#### Multicolinearidade ####

# Atribuir semente
set.seed(1)

# Criar dados simulados de área, número de cômodos e preço dos imóveis
area <- rnorm(100, mean = 1500, sd = 500)
comodos <- round(rnorm(100, mean = 3 + 0.5 * area/100, sd = 1), 0)
preco <- rnorm(100, mean = 200000 + 100 * area + 50000 * comodos, sd = 50000)
area <- area / 10

# Juntar os dados em um data frame
dados <- data.frame(area, comodos, preco)

# Criar gráfico que demonstra multicolinearidade entre área e cômodos
ggplot(dados, aes(x = area, y = preco, color = comodos)) +
  geom_point() +
  scale_color_viridis_c(name = "Cômodos") +
  scale_y_continuous(labels = scales::dollar_format(prefix = "R$ ",
                                                    big.mark = ".")) +
  scale_x_continuous(labels = \(x) unname(latex2exp::TeX(paste0(
    as.character(x), "$m^2$"
  )))) +
  theme_minimal(20) +
  labs(x = "Área",
       y = "Preço do imóvel")


# Carregar pacote `car` para calcular variance inflation factor (VIF)
library(car)

# Ajustar modelo do preço do imóvel explicado pela área e número de cômodos
regressao_linear <- lm(preco ~ area + comodos,
                       data = dados)

# Calcular VIF dos preditores do modelo
vif(regressao_linear)

#### Outliers ####

# Atribuir semente
set.seed(1)

# Criar dados simulados e juntar em um data frame
x <- 1:10
y <- 2*x + rnorm(10, mean = 0, sd = 1)
df <- data.frame(x, y)

# Criando o gráfico sem outlier
p_sem_outlier <- ggplot(df, aes(x = x, y = y)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  scale_x_continuous(limits = c(0, 12.5)) +
  scale_y_continuous(limits = c(0, 50)) +
  theme_minimal(18) +
  labs(x = unname(latex2exp::TeX("$x_1$")),
       y = unname(latex2exp::TeX("$Y$")))

# Adicionando um outlier aos dados
df[11, ] <- c(12, 50)

# Criando o gráfico com outlier
p_com_outlier <- ggplot(df, aes(x = x, y = y)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  theme_minimal(18) +
  scale_x_continuous(limits = c(0, 12.5)) +
  scale_y_continuous(limits = c(0, 50)) +
  labs(x = unname(latex2exp::TeX("$x_1$ com \\textit{outlier}")),
       y = NULL)

# Juntar ambos os gráficos com `patchwork`
p_sem_outlier + p_com_outlier

# Criar modelo com dados que possuem outlier
modelo_com_outlier <- lm(y ~ x, data = df)

# Calcular a Distância de Cook para cada observação
cooks.distance(modelo_com_outlier)

# Calcular a medida de 4/n para comparar os valores de Distância de Cook calculados acima
4 / nrow(df)

# Atribuir as Distâncias de Cook em um objeto
distancia_cook <- cooks.distance(modelo_com_outlier)

# Criar data frame com as Distâncias para cada observação
df_cook <- data.frame(distancia_cook, row.names = NULL)

# Criar gráfico que mostra como comparar a Distância de Cook ao ponto de corte (4/n)
ggplot(df_cook, aes(x = as.numeric(row.names(df_cook)), y = distancia_cook)) +
  geom_point(size = 3) +
  scale_x_continuous(breaks = 1:11) +
  geom_hline(
    yintercept = 4 / nrow(df_cook),
    color = "darkred",
    size = 2
  ) +
  labs(x = "Observação", y = "Distância de Cook") +
  theme_minimal(20)

```

## Regressõa logistica binomial/multinomial
```{r}
# Regressao logistica multinomial 
# Baseado em https://fileds.quarto.pub/data-analytics-with-r/logistic_regression_multiple.html
# Dados de pacientes com diabetes para avaliacao de preditores

library(tidyverse) 

# install.packages("mlbench") # Install to get PimaIndiansDiabetes2 dataset
library(mlbench)
data("PimaIndiansDiabetes2") # Using 2 because it is a corrected version.
View(PimaIndiansDiabetes2)
print(PimaIndiansDiabetes2)

## Visualizacao de graficos pareados
library(patchwork)

p1 <- PimaIndiansDiabetes2 |>
  ggplot(aes(y = pressure)) + 
  facet_grid(~ diabetes) +
  geom_boxplot() +
  theme(
    axis.text.x = element_blank(),  # Remove axis text
    axis.ticks.x = element_blank()  # Remove axis title
  ) +
  labs(y = "Diastolic blood pressure (mmHg)")

p2 <- PimaIndiansDiabetes2 |>
  ggplot(aes(x = pressure, fill = diabetes)) + 
  geom_density(alpha = 0.5) +
  theme(
    legend.position = "top",
    axis.text.y = element_blank(),  # Remove axis text
    axis.ticks.y = element_blank()  # Remove axis title
  ) +
  labs(
    x = "Diastolic blood pressure (mmHg)",
    y = "",
    fill = "Diabetes"
  )

# patchwork para pareamento 
p1 + p2

## Parece existir uma relacao entre o aumento da pressao sanguinea e a diabete. Mas nao e claro

## Modelo logistico binomial
model <- glm(diabetes ~ pressure, data = PimaIndiansDiabetes2, family = "binomial")
model
model_summary <- model |> summary()
model_summary$coefficients

# Relação positiva e estatisticamente significativa entre pressão arterial e o fato de ter diabetes
# ao nível de significância de 0,05.
# E se incluirmos mais variáveis em nosso modelo, poderia haver alguma outra explicação para a 
# relação entre pressão arterial e diabetes?

# Avaliar a multicolinearidade
# Instalar se necessário
# install.packages("corrplot")

library(corrplot)

# Matriz de correlação das variáveis numéricas 
#sempre avaliar a presenca de NAs
corr_matrix <- cor(PimaIndiansDiabetes2[,1:8])

## Remove os valores faltantes
corr_matrix <- cor(PimaIndiansDiabetes2[,1:8], use = "complete.obs")

# Gráfico da correlação
corrplot(corr_matrix, method = "color", 
         type = "upper",      # apenas a metade superior
         tl.col = "black",    # cor dos rótulos
         addCoef.col = "white", # adiciona os valores da correlação
         number.cex = 0.7, 
         tl.srt = 45)         # rotação dos rótulos

### Visualizacao 2 da correlacao
# install.packages("ggcorrplot")
library(ggcorrplot)
?ggcorrplot
ggcorrplot(corr_matrix, 
           hc.order = TRUE,  # ordena hierarquicamente
           type = "upper",
           lab = TRUE,       # mostra os números
           lab_size = 3, 
           colors = c("blue", "white", "red"))

# Valores próximos de ±1 indicam alta correlação.
# Mas a correlação simples só captura relações lineares entre pares de variáveis, não multivariadas.

##### Avaliacao da Variance Inflation Factor (VIF)

####### Primeiro gerar o modelo 
model1 <- PimaIndiansDiabetes2 |> 
  glm(formula = diabetes ~ pressure + mass, family = "binomial")
model1 |> summary()

# Instalar se necessário
# install.packages("car")

library(car)

vif(model1)

# O VIF é o método mais usado para detectar multicolinearidade. Valores comuns de referência:
#   
#   VIF > 5 → possível multicolinearidade
# 
#   VIF > 10 → multicolinearidade alta

#### Interpretar o modelo 

summary(model1)
# Observamos que tanto o índice de massa corporal (mass)
# quanto a pressão arterial (pressure) são estatisticamente significativos
# quando incluídos no modelo e que ambos aumentam a probabilidade de diabetes,
# ou seja, seus coeficientes são positivos

#######
"Também podemos optar por incluir a variável triceps, que mede a espessura da prega 
cutânea do tríceps e que se correlaciona com a gordura corporal. Vamos ajustar
um modelo usando triceps em vez de mass."
####

model2 <- PimaIndiansDiabetes2 |> 
  glm(formula = diabetes ~ pressure + triceps, family = "binomial")
model2 |> summary()

#####
"Ambas as variáveis são significativas ao nível de 0,05. 
Mas qual modelo é o melhor?
Na regressão linear usamos o R² 
Para regressão logística usamos o Critério de Informação de Akaike (AIC).
O AIC depende do número de observações, então, diferentemente do R², 
o valor isolado do AIC de um modelo não é muito informativo. Ele é útil para comparar
modelos e decidir qual modelo é mais adequado.
No resumo do modelo (summary(model)), a penúltima linha mostra o AIC. 
Ele estima o erro de predição do modelo, sendo que valores menores de AIC indicam modelos melhores."

# para acessar rapidamente diretamente no modelo
model1$aic
model2$aic

## Lembrar de ver a VIF
vif(model2)

### Exercicio para avaliar depois 
#fazer um modelo completo com os tres dados avaliados e compara o IAC mass triceps e pressure
names(PimaIndiansDiabetes2)

model3 <- PimaIndiansDiabetes2 |> 
  glm(formula = diabetes ~ pressure + triceps + mass, family = "binomial")
model3 |> summary()

# para acessar rapidamente 
model1$aic
model2$aic
model3$aic

#### Selecao de modelo
# Em regressão logística múltipla, podemos usar diferentes estratégias:
# - Conhecimento de domínio (D/E knowledge): escolher variáveis baseadas em
#   experiência clínica ou científica.
# - Métodos automáticos:
#    * backward: começa com todas as variáveis, remove as menos significativas
#    * forward: começa com nenhuma variável, adiciona as mais significativas
#    * stepwise/mista: combina adição e remoção de variáveis
# A escolha do modelo final pode ser baseada em:
# - Significância das variáveis (p < 0.05)
# - Critério de Informação de Akaike (AIC)
# - Conhecimento de domínio
#Podemos incluir todas as variáveis do dataset usando o ponto .
#e construir um modelo usando seleção backward:

# Usando "." para incluir todas as variáveis do dataset
model_full <- glm(diabetes ~ ., data = PimaIndiansDiabetes2, family = "binomial")

summary(model_full)

"temos varias variaives que nao foram significativas. A pressao foi a que teve o maior p
(p>0.9) podemos criar um modelo sem ela para avaliarmos" 

model <- PimaIndiansDiabetes2 |> 
  glm(formula = diabetes ~ . - pressure, family = "binomial")
model_summary <- model |> summary()
model_summary$aic

# Fazendo de forma automatica

"lembrar qeu para isso, e necessario remover os NAs"
df_complete <- na.omit(PimaIndiansDiabetes2)

model_full <- glm(diabetes ~ ., data = df_complete, family = "binomial")
summary(model_full)
# Seleção backward
model_selected <- step(model_full, direction = "backward")
model_selected <- step(model_full, direction = "forward")
model_selected <- step(model_full, direction = "both")

# O step() remove variáveis de forma iterativa para minimizar o AIC
model_selected <- step(model_full, direction = "backward")
#Analisando o modelo final
summary(model_selected)

# e se tirarmos a gravidez do modelo, sera que nao ficaria melhor ?????
glm(formula = diabetes ~ glucose + mass + pedigree + 
      age, family = "binomial", data = df_complete)

### Predicao 

# Assim como na regressão logística simples, podemos usar nosso modelo
# para **prever probabilidades**. Para isso, precisamos fornecer um
# data.frame (ou tibble) com valores para os preditores incluídos no modelo.

# Exemplo: modelo com triceps e pressão arterial
model_triceps <- PimaIndiansDiabetes2 |> 
  glm(formula = diabetes ~ pressure + triceps, family = "binomial")

# Queremos prever a probabilidade de diabetes em 5 anos para uma mulher
# com pressão diastólica = 75 mm Hg e espessura da prega cutânea do tríceps = 30 mm
prediction_prob <- model_triceps |> 
  predict(data.frame(pressure = 75, triceps = 30), type = "response")

# O parâmetro type = "response" retorna a **probabilidade** (entre 0 e 1)
prediction_prob
# Exemplo de saída: 0.3471786
# Interpretação: segundo nosso modelo, a probabilidade de desenvolver diabetes
# nesta situação é aproximadamente 34,7%.

# ---------------------------------------------------------------
# Transformando probabilidade em classe
# ---------------------------------------------------------------
# Se quisermos classificar como "diabetes positivo" ou "diabetes negativo",
# podemos usar um **cutoff de 0.5**:
# - probabilidade > 0.5 → "positivo"
# - probabilidade <= 0.5 → "negativo"

prediction_class <- ifelse(prediction_prob > 0.5, "pos", "neg")
prediction_class
# Saída: "neg"
# Interpretação: neste caso, a previsão é negativa, pois a probabilidade
# ficou abaixo de 0.5.


# ---------------------------------------------------------------
# Avaliando o Desempenho do Modelo – Accuracy
# ---------------------------------------------------------------

# Além do AIC, podemos avaliar o desempenho do modelo usando **predição**,
# calculando a **accuracy** (acurácia) do modelo.
# Accuracy = número de predições corretas / número total de predições

# Vamos calcular a acurácia para todas as observações completas do dataset
model <- PimaIndiansDiabetes2 |> 
  glm(formula = diabetes ~ . - pressure, family = "binomial")

# Passo a passo:
accuracy_result <- PimaIndiansDiabetes2 |> 
  # 1️⃣ Calcular a probabilidade prevista de diabetes para todas as observações
  mutate(predicted_prob = model |> predict(PimaIndiansDiabetes2, type = "response")) |> 
  # 2️⃣ Transformar probabilidade em classe: >0.5 = "pos", <=0.5 = "neg"
  mutate(predicted_outcome = ifelse(predicted_prob > 0.5, "pos", "neg")) |> 
  # 3️⃣ Remover casos com valores ausentes (se houver)
  filter(!is.na(predicted_outcome)) |> 
  # 4️⃣ Calcular a acurácia: proporção de predições corretas
  summarize(accuracy = mean(predicted_outcome == diabetes))

accuracy_result
# Exemplo de saída: 0.7831633
# Interpretação: o modelo acerta aproximadamente 78% das observações.

# ---------------------------------------------------------------
# Observações sobre a acurácia
# ---------------------------------------------------------------
# - Accuracy é intuitiva, mas pode ser enganosa se as classes forem desbalanceadas.
# - Por exemplo, se o modelo prevê melhor os casos negativos do que positivos,
#   a acurácia global não mostra essa diferença.
# - Podemos visualizar a performance por classe (diabetes positivo e negativo)
#   para entender melhor onde o modelo acerta e erra.

########## Visualizacao da acuracia 

library(ggplot2)
library(dplyr)
library(scales)
set.seed(42)

PimaIndiansDiabetes2 %>%
  mutate(predicted_prob = model |> predict(PimaIndiansDiabetes2, type = "response")) %>%
  drop_na() %>%
  mutate(predicted_outcome = factor(ifelse(predicted_prob > 0.5, "pos", "neg"))) %>%
  mutate(correct_prediction = factor(ifelse(predicted_outcome == diabetes, "yes", "no"))) %>%
  mutate(index = row_number(), .by = diabetes) %>%
  arrange(index) %>%
  ggplot(aes(x = diabetes, y = index, col = correct_prediction)) +
  geom_jitter(width = 0.2, height = 0, size = 2, alpha = 0.8) +
  scale_color_manual(values = c("yes" = "#1b9e77", "no" = "#d95f02")) +
  labs(
    title = "Visualização da Acurácia do Modelo",
    subtitle = "Cada ponto representa uma observação do dataset",
    x = "Classe observada (Diabetes)",
    y = "Observações (index)",
    col = "Predição correta"
  ) +
  theme_bw() +
  theme(
    text = element_text(size = 14),   # tamanho da fonte
    legend.position = "top",
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )


#########
# install.packages("pROC")
library(pROC)

# Selecionar linhas completas para as variáveis do modelo
complete_rows <- complete.cases(PimaIndiansDiabetes2[, c("diabetes", "pressure", "triceps")])
data_complete <- PimaIndiansDiabetes2[complete_rows, ]

# Ajustar o modelo apenas com dados completos
model <- glm(diabetes ~ pressure + triceps, data = data_complete, family = "binomial")

# Obter probabilidades previstas
pred_prob <- predict(model, type = "response")

# Curva ROC
roc_curve <- roc(data_complete$diabetes, pred_prob)

# Plotar
plot(roc_curve, col = "#1b9e77", lwd = 2, main = "Curva ROC - Modelo de Diabetes")
auc(roc_curve)  # Área sob a curva
###
"Quanto mais a curva verde se aproxima do canto superior esquerdo, 
melhor o modelo separa os casos positivos dos negativos.
Se a curva segue próxima da diagonal, o modelo não tem habilidade de discriminação.
A curva está acima da diagonal, indicando que o modelo 
tem capacidade moderada de separar classes."

## interpretaçao
"AUC = 0.5 → o modelo não consegue separar as classes melhor que o acaso.
AUC = 1 → o modelo separa perfeitamente as classes.
AUC = 0.6768 → o modelo tem capacidade moderada de discriminar entre positivos e negativos.
Em termos simples: se pegássemos aleatoriamente uma pessoa com diabetes e uma sem diabetes,
nosso modelo acertaria quem tem diabetes aproximadamente 67,7% das vezes."

```

## Regressão polinomial
```{r}
# =========================================
# REGRESSÃO POLINOMIAL - SCRIPT DIDÁTICO
# =========================================
# Objetivo: Ajustar polinômios de diferentes graus,
# comparar os modelos e identificar o ponto ótimo
# =========================================

# -----------------------------
# 0. Carregar pacotes necessários
# -----------------------------
# install.packages("tidyverse")  # Execute se ainda não tiver
library(tidyverse)

# -----------------------------
# 1. Gerar dados simulados
# -----------------------------

set.seed(123)             
n <- 50                   
dose <- seq(0, 10, length.out = n)  

# Resposta simulada: efeito quadrático descendente + ruído
# Começa em torno de 120 mmHg, atinge um pico, depois decresce
pressure <- 120 - 5*dose + 0.4*dose^2 + rnorm(n, sd=2)  

data <- data.frame(dose, pressure)

# Visualizar os dados
ggplot(data, aes(x=dose, y=pressure)) +
  geom_point(color="darkgreen", size=2) +
  labs(title="Efeito da Dose sobre Pressão Arterial",
       x="Dose (unidades)", y="Pressão (mmHg)") +
  theme_classic()

# Visualizar os dados com linha de regressão linear
ggplot(data, aes(x=dose, y=pressure)) +
  geom_point(color="darkgreen", size=2) +                  # Pontos observados
  geom_smooth(method="lm", formula = y ~ x, color="blue", se=FALSE, size=1) +  # Regressão linear
  labs(title="Efeito da Dose sobre Pressão Arterial",
       x="Dose (unidades)", y="Pressão (mmHg)") +
  theme_classic()

# -----------------------------
# 2. Ajustar polinômios de diferentes graus
# -----------------------------
# Modelo Linear (grau 1)
model1 <- lm(pressure ~ dose, data=data)

# Modelo Quadrático (grau 2)
model2 <- lm(pressure ~ poly(dose, 2, raw=TRUE), data=data)

# Modelo Cúbico (grau 3)
model3 <- lm(pressure ~ poly(dose, 3, raw=TRUE), data=data)

# -----------------------------
# 3. Comparar os modelos
# -----------------------------
# Avaliar R² (quanto do efeito da dose é explicado)
r2_values <- c(
  Linear = summary(model1)$r.squared,
  Quadratico = summary(model2)$r.squared,
  Cubico = summary(model3)$r.squared
)
print("R² de cada modelo:")
print(r2_values)

# Comparar AIC (quanto penaliza complexidade do modelo)
print("AIC dos modelos:")
print(AIC(model1, model2, model3))

# Observações:
# - R² mais alto → modelo explica melhor a variação da resposta
# - AIC mais baixo → melhor ajuste penalizando complexidade

# -----------------------------
# 4. Visualizar ajuste dos polinômios com legenda automática
# -----------------------------

# Criar um dataframe "long" com predições de cada modelo
data_pred <- data %>%
  mutate(
    Linear = predict(model1, newdata = data),
    Quadratico = predict(model2, newdata = data),
    Cubico = predict(model3, newdata = data)
  ) %>%
  pivot_longer(cols = c("Linear", "Quadratico", "Cubico"),
               names_to = "Modelo", values_to = "Predicao")

# Plotando
ggplot(data, aes(x=dose, y=pressure)) +
  geom_point(size=2, color="black") +  # pontos observados
  geom_line(data=data_pred, aes(x=dose, y=Predicao, color=Modelo), size=1.2) +
  scale_color_manual(values=c("Linear"="dodgerblue", 
                              "Quadratico"="firebrick", 
                              "Cubico"="darkgreen")) +
  labs(title="Comparação de Modelos Polinomiais",
       x="Dose (unidades)", y="Pressão (mmHg)", color="Modelo") +
  theme_classic() +
  theme(legend.position = "top", legend.title = element_text(face="bold"))

# -----------------------------
# 5. Encontrar ponto de máximo/minimo (polinômio quadrático)
# -----------------------------
# Fórmula do polinômio: y = a + b*x + c*x^2
coefs <- coef(model2)   # coeficientes: Intercept, dose, dose^2
a <- coefs[1]
b <- coefs[2]
c <- coefs[3]
"Aqui estamos trabalhando com um polinômio quadrático, que tem a forma geral:
y = a + b*x + c*x^2
y → variável resposta (pressão arterial).
x → variável preditora (dose do medicamento).
a → intercepto, valor de y quando x=0.
b → coeficiente linear, indica a inclinação inicial da curva.
c → coeficiente quadrático, determina a curvatura (para cima ou para baixo).
No caso de um polinômio quadrático, a curva é uma parábola que pode ter um máximo 
(curva voltada para baixo, c < 0) ou mínimo (curva voltada para cima, c > 0)."

# Ponto ótimo (x*) = -b / (2*c)
x_opt <- -b / (2*c)

# Valor da resposta neste ponto
y_opt <- a + b*x_opt + c*x_opt^2

cat("Dose ótima:", round(x_opt,2), "\n")
cat("Pressão prevista no ponto ótimo:", round(y_opt,2), "\n")

# -----------------------------
# 6. Visualizando ponto ótimo no gráfico
# -----------------------------
ggplot(data, aes(x=dose, y=pressure)) +
  geom_point(size=2) +
  geom_smooth(method="lm", formula = y ~ poly(x,2,raw=TRUE), color="red", se=FALSE) +
  geom_point(aes(x=x_opt, y=y_opt), color="blue", size=4) +
  geom_text(aes(x=x_opt, y=y_opt, label=paste0("Dose óptima = ", round(x_opt,2))), 
            vjust=-1, hjust=0.5, color="blue", size=5) +
  labs(title="Ponto Ótimo do Polinômio Quadrático",
       x="Dose (unidades)", y="Pressão (mmHg)") +
  theme_minimal()

# -----------------------------
# 6. Encontrar pontos críticos (cúbico)
# -----------------------------
coefs3 <- coef(model3)  # Intercept, dose, dose^2, dose^3
a <- coefs3[1]; b <- coefs3[2]; c <- coefs3[3]; d <- coefs3[4]

# Derivada do cúbico: dy/dx = b + 2*c*x + 3*d*x^2 = 0

# Resolver a derivada usando polyroot
raizes <- polyroot(c(b, 2*c, 3*d))   # retorna todas as raízes (complexas)
# Filtrar apenas as reais
pontos_criticos <- Re(raizes[Im(raizes) == 0])

# Filtrar apenas valores dentro do intervalo da dose
pontos_criticos <- pontos_criticos[pontos_criticos >= min(dose) & pontos_criticos <= max(dose)]

# Calcular y nos pontos críticos
y_criticos <- a + b*pontos_criticos + c*pontos_criticos^2 + d*pontos_criticos^3

# Mostrar resultados (somente se houver pontos críticos)
if(length(pontos_criticos) > 0){
  df_criticos <- data.frame(
    dose = pontos_criticos,
    pressure = y_criticos
  )
  cat("\nPolinômio cúbico - pontos críticos dentro do intervalo da dose:\n")
  print(df_criticos)
} else {
  df_criticos <- NULL
  cat("\nNão há pontos críticos reais dentro do intervalo da dose.\n")
}

# -----------------------------
# 7. Visualizar pontos críticos no gráfico (cúbico)
# -----------------------------
ggplot(data, aes(x=dose, y=pressure)) +
  geom_point(size=2, color="black") +  # pontos observados
  geom_line(data=data_pred, aes(x=dose, y=Predicao, color=Modelo), size=1.2) +
  # Adicionar pontos críticos somente se existirem
  {if(!is.null(df_criticos)) geom_point(data=df_criticos, aes(x=dose, y=pressure), color="blue", size=4)} +
  {if(!is.null(df_criticos)) geom_text(data=df_criticos, aes(x=dose, y=pressure, 
                                                             label=paste0("x=",round(dose,2))),
                                       vjust=-1, color="blue", size=4)} +
  scale_color_manual(values=c("Linear"="dodgerblue", 
                              "Quadratico"="firebrick", 
                              "Cubico"="darkgreen")) +
  labs(title="Polinômio Cúbico: pontos críticos",
       x="Dose (unidades)", y="Pressão (mmHg)", color="Modelo") +
  theme_classic() +
  theme(legend.position = "top", legend.title = element_text(face="bold"))


####
"Nem sempre o modelo cúbico vai gerar pontos críticos visíveis 
nos dados que você simulou.
Nem todo modelo mais complexo garante “picos” dentro do intervalo analisado.
Mesmo sem pontos críticos, o gráfico da curva cúbica ainda é útil para comparar o ajuste com os modelos linear e quadrático."


```