# 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
## qual teste estatístico usar
# 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
Ir em instalar pela primeira vez
Escolher a versão mais atual
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
# 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**
## 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.
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."
```
Sign in with Wallet
Connect another wallet