Regressão Linear Múltipla — Análise Completa da Base Iris

Estatística Aplicada — Universidade Federal de Campina Grande (UFCG)

Ideia central. Regressão múltipla não é uma sequência de botões. É um argumento científico: definimos uma pergunta, explicitamos um modelo, confrontamos suas suposições com os dados, avaliamos incerteza e só então comunicamos o que o modelo permite — e o que ele não permite — concluir.

1 Visão geral da análise

1.1 Preparação do ambiente

Instale os pacotes uma única vez. Depois, abra o arquivo .Rmd no RStudio e use Knit, ou execute rmarkdown::render(). No HTML, os códigos ficam recolhidos: clique em Code para visualizá-los. Para reproduzir ou alterar a análise, execute/renderize o documento inteiro na ordem apresentada, pois cada etapa usa objetos criados nas etapas anteriores.

install.packages(c(
  "rmarkdown", "knitr", "ggplot2", "dplyr", "tidyr", "purrr",
  "broom", "performance", "patchwork", "GGally",
  "kableExtra", "DT"
))

rmarkdown::render("tutorial_regressao_multipla_iris.Rmd")

1.2 Objetivos da análise

Esta análise foi estruturada para:

1. transformar uma pergunta substantiva em uma equação de regressão;
2. distinguir variável resposta, preditores numéricos e preditores categóricos;
3. realizar uma análise exploratória orientada ao modelo;
4. ajustar e interpretar uma regressão linear múltipla em R;
5. verificar linearidade, independência, homocedasticidade, normalidade aproximada, multicolinearidade e influência;
6. comparar cinco modelos usando teste F, \(R^2\) ajustado, AIC, BIC e parcimônia;
7. interpretar corretamente variáveis categóricas e categorias de referência;
8. diferenciar intervalo de confiança da média e intervalo de previsão individual;
9. reconhecer multicolinearidade, extrapolação e linguagem causal indevida;
10. selecionar um modelo com base no conjunto das evidências e produzir uma conclusão defensável.

1.3 O mapa do processo

Etapa	Pergunta que guia a etapa	Produto
1. Problema	O que queremos explicar ou prever?	pergunta, população e unidade de análise
2. Dados	As variáveis medem o que precisamos?	dicionário, qualidade e exploração
3. Especificação	Que relação estamos supondo?	fórmula, termos e interações justificadas
4. Ajuste	Quais coeficientes melhor reproduzem os dados?	modelo estimado e incerteza
5. Comparação preliminar	Qual modelo tem melhor ajuste e quais termos acrescentam informação?	tabela comparativa, teste F e candidato ao diagnóstico
6. Diagnóstico	O candidato respeita os pressupostos de forma razoável?	resíduos, influência, VIF e gráficos auxiliares
7. Decisão	Qual especificação é mais defensável para cada objetivo?	melhor ajuste, alternativa interpretável e justificativa
8. Comunicação	O que aprendemos e com quais limites?	conclusão, intervalos e próximos passos

2 A pergunta vem antes do modelo

Usaremos Sepal.Length — o comprimento da sépala, em centímetros — como resposta. Queremos investigar e prever essa medida a partir da largura da sépala, das medidas da pétala e da espécie.

Uma pergunta bem formulada seria:

Entre flores semelhantes nas demais características consideradas, como o comprimento esperado da sépala varia com cada preditor, e quão bem o modelo prevê novas flores?

Essa frase já contém duas metas diferentes:

• inferência/explicação: interpretar relações condicionais e sua incerteza;
• previsão: estimar valores para novas observações com baixo erro.

“Melhor modelo” não é uma propriedade absoluta. Um modelo pode ser excelente para previsão e ruim para interpretação; outro pode ser transparente e estável, mas perder um pouco de acurácia. Antes de comparar modelos, declare o objetivo.

2.1 O modelo matemático

Para a observação \(i\), um modelo com três preditores numéricos é

\[ Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \beta_3 X_{3i} + \varepsilon_i, \]

em que:

• \(Y_i\) é o comprimento observado da sépala;
• \(\beta_0\) é o valor esperado de \(Y\) quando os preditores valem zero;
• \(\beta_j\) é a mudança média esperada em \(Y\) para uma unidade adicional de \(X_j\), mantendo os demais preditores constantes;
• \(\varepsilon_i\) representa tudo o que o modelo não explica.

Quando incluímos Species, o R cria contrastes (as chamadas variáveis indicadoras ou dummies) comparando cada nível com uma categoria de referência.

2.2 Suposições: o que realmente precisamos verificar

Suposição	Significado operacional	Como verificar	O que fazer se falhar
Linearidade	a média condicional de \(Y\) é linear nos parâmetros	resíduos vs. ajustados; resíduos parciais	transformação, polinômio, spline ou interação justificada
Independência	erros de unidades distintas não são correlacionados	desenho da coleta; ordem/tempo/grupos	modelos mistos, séries temporais, erros agrupados
Variância aproximadamente constante	dispersão dos erros é semelhante ao longo da média	resíduos vs. ajustados; escala-localização	transformação, erros robustos, WLS
Normalidade aproximada dos erros	necessária principalmente para inferência exata em amostras pequenas	QQ-plot; forma dos resíduos	bootstrap, erros robustos, investigar modelo/outliers
Ausência de colinearidade prejudicial	preditores não carregam quase a mesma informação	VIF/GVIF; estabilidade dos coeficientes	reformular, combinar, centrar, regularizar, coletar mais dados
Ausência de influência excessiva	nenhuma observação domina o ajuste	alavancagem e distância de Cook	verificar dado, relatar sensibilidade; nunca excluir automaticamente
Especificação adequada	termos importantes não foram omitidos	teoria, resíduos, interações, validação	reespecificar com justificativa substantiva

Cheque conceitual. Normalidade dos preditores não é pressuposto da regressão linear. O foco é a distribuição dos erros condicionais, e nem toda pequena violação invalida o modelo.

3 Conhecendo a base Iris

3.1 Estrutura, unidade e qualidade

dados <- datasets::iris |>
  mutate(
    Species = factor(Species),
    id = row_number()
  )

glimpse(dados)

## Rows: 150
## Columns: 6
## $ Sepal.Length <dbl> 5.1, 4.9, 4.7, 4.6, 5.0, 5.4, 4.6, 5.0, 4.4, 4.9, 5.4, 4.…
## $ Sepal.Width  <dbl> 3.5, 3.0, 3.2, 3.1, 3.6, 3.9, 3.4, 3.4, 2.9, 3.1, 3.7, 3.…
## $ Petal.Length <dbl> 1.4, 1.4, 1.3, 1.5, 1.4, 1.7, 1.4, 1.5, 1.4, 1.5, 1.5, 1.…
## $ Petal.Width  <dbl> 0.2, 0.2, 0.2, 0.2, 0.2, 0.4, 0.3, 0.2, 0.2, 0.1, 0.2, 0.…
## $ Species      <fct> setosa, setosa, setosa, setosa, setosa, setosa, setosa, s…
## $ id           <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17…

Cada linha representa uma flor; cada medida está em centímetros. Há 150 observações, 50 de cada espécie.

dicionario <- tibble::tribble(
  ~variavel, ~tipo, ~papel, ~definicao,
  "Sepal.Length", "numérica", "resposta", "comprimento da sépala (cm)",
  "Sepal.Width", "numérica", "preditor", "largura da sépala (cm)",
  "Petal.Length", "numérica", "preditor", "comprimento da pétala (cm)",
  "Petal.Width", "numérica", "preditor", "largura da pétala (cm)",
  "Species", "categórica", "preditor/estrato", "espécie da flor",
  "id", "inteira", "identificador", "linha original para rastreabilidade"
)

knitr::kable(dicionario, align = c("l", "l", "l", "l")) |>
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

variavel	tipo	papel	definicao
Sepal.Length	numérica	resposta	comprimento da sépala (cm)
Sepal.Width	numérica	preditor	largura da sépala (cm)
Petal.Length	numérica	preditor	comprimento da pétala (cm)
Petal.Width	numérica	preditor	largura da pétala (cm)
Species	categórica	preditor/estrato	espécie da flor
id	inteira	identificador	linha original para rastreabilidade

qualidade <- dados |>
  summarise(
    observacoes = n(),
    duplicatas_exatas = sum(duplicated(across(-id))),
    valores_ausentes = sum(is.na(across(everything()))),
    especies = n_distinct(Species)
  )

qualidade

Há uma duplicata exata nas medidas e espécie. Isso não prova erro: duas flores podem ter medidas idênticas. Sem um identificador original da coleta, não devemos apagar a linha automaticamente. O correto é registrar a ambiguidade e fazer análise de sensibilidade se ela puder afetar o resultado.

3.2 Estatísticas descritivas

descritivas <- dados |>
  select(-id) |>
  pivot_longer(where(is.numeric), names_to = "variavel", values_to = "valor") |>
  group_by(variavel) |>
  summarise(
    n = n(),
    media = mean(valor),
    dp = sd(valor),
    minimo = min(valor),
    q1 = quantile(valor, 0.25),
    mediana = median(valor),
    q3 = quantile(valor, 0.75),
    maximo = max(valor),
    .groups = "drop"
  )

knitr::kable(descritivas, digits = 2) |>
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

variavel	n	media	dp	minimo	q1	mediana	q3	maximo
Petal.Length	150	3.76	1.77	1.0	1.6	4.35	5.1	6.9
Petal.Width	150	1.20	0.76	0.1	0.3	1.30	1.8	2.5
Sepal.Length	150	5.84	0.83	4.3	5.1	5.80	6.4	7.9
Sepal.Width	150	3.06	0.44	2.0	2.8	3.00	3.3	4.4

3.3 Distribuições e grupos

dados_longos <- dados |>
  pivot_longer(
    cols = c(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width),
    names_to = "medida",
    values_to = "valor"
  )

ggplot(dados_longos, aes(Species, valor, color = Species)) +
  geom_boxplot(aes(fill = Species), alpha = 0.12, width = 0.58,
               outlier.shape = NA, linewidth = 0.55) +
  geom_jitter(width = 0.14, alpha = 0.55, size = 1.35) +
  facet_wrap(~ medida, scales = "free_y", ncol = 2) +
  scale_color_manual(values = paleta_especies) +
  scale_fill_manual(values = paleta_especies) +
  labs(
    title = "Medidas morfológicas por espécie",
    subtitle = "150 flores; 50 observações por espécie; valores em centímetros",
    x = NULL, y = "Valor (cm)", color = "Espécie", fill = "Espécie"
  ) +
  guides(fill = "none")

Distribuições das quatro medidas por espécie. Pontos e boxplots tornam a amostra e a dispersão visíveis sem depender apenas das médias.

Leitura. A espécie separa fortemente as medidas da pétala e também parte do comprimento da sépala. Portanto, uma relação agregada entre pétala e sépala pode misturar dois fenômenos: diferenças entre espécies e variação dentro de cada espécie.

3.4 Relações bivariadas não contam a história completa

GGally::ggpairs(
  dados,
  columns = 1:4,
  mapping = aes(color = Species, alpha = 0.7),
  upper = list(continuous = GGally::wrap("cor", size = 3.4)),
  diag = list(continuous = GGally::wrap("densityDiag", alpha = 0.35))
) +
  scale_color_manual(values = paleta_especies) +
  scale_fill_manual(values = paleta_especies) +
  theme_bw(base_size = 10)

Matriz de pares: distribuições na diagonal, correlações no triângulo superior e dispersões no triângulo inferior.

matriz_cor <- cor(select(dados, Sepal.Length:Petal.Width))
round(matriz_cor, 3)

##              Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length        1.000      -0.118        0.872       0.818
## Sepal.Width        -0.118       1.000       -0.428      -0.366
## Petal.Length        0.872      -0.428        1.000       0.963
## Petal.Width         0.818      -0.366        0.963       1.000

Petal.Length e Petal.Width têm correlação muito alta. Isso sinaliza possível multicolinearidade, mas a decisão não deve ser tomada por um limiar mecânico de correlação. O VIF, a estabilidade dos coeficientes e o objetivo do modelo completam o diagnóstico.

Correção de uma armadilha comum. Dois preditores positivamente correlacionados não precisam ter coeficientes múltiplos com o mesmo sinal. Cada coeficiente é um efeito parcial, condicionado aos demais termos. Uma inversão de sinal pode revelar colinearidade, confundimento ou uma pergunta condicional diferente — não uma contradição lógica automática.

p_agregado <- ggplot(dados, aes(Petal.Length, Sepal.Length)) +
  geom_point(color = azul, alpha = 0.65) +
  geom_smooth(method = "lm", se = TRUE, color = azul_escuro, fill = "#DCE8F7") +
  labs(
    title = "Relação agregada",
    subtitle = "Uma única reta para todas as flores",
    x = "Comprimento da pétala (cm)", y = "Comprimento da sépala (cm)"
  )

p_grupos <- ggplot(dados, aes(Petal.Length, Sepal.Length, color = Species)) +
  geom_point(alpha = 0.72) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 0.9) +
  scale_color_manual(values = paleta_especies) +
  labs(
    title = "Relações por espécie",
    subtitle = "A inclinação pode variar entre grupos",
    x = "Comprimento da pétala (cm)", y = "Comprimento da sépala (cm)",
    color = "Espécie"
  )

p_agregado + p_grupos

A linha agregada mistura diferenças entre espécies e relações dentro de cada espécie.

4 Construção progressiva dos cinco modelos

Os cinco modelos serão ajustados com as 150 observações. A ordem foi escolhida para que cada novo ajuste responda a uma dúvida surgida no modelo anterior. A decisão será tomada somente depois que todos tiverem sido estimados, interpretados e comparados.

4.1 Como ler a saída de lm()

• Estimate: estimativa de cada coeficiente \(\beta_j\);
• Std. Error: erro-padrão, que mede a incerteza da estimativa;
• t value e p-valor: teste de \(H_0: \beta_j=0\), condicionado aos demais termos do modelo;
• Residual standard error: escala típica dos resíduos, em centímetros;
• Multiple \(R^2\): proporção da variabilidade observada explicada pelo ajuste;
• Adjusted \(R^2\): versão do \(R^2\) que considera o número de parâmetros;
• F-statistic: teste global do modelo contra um ajuste somente com intercepto.

O coeficiente estimado e seu erro-padrão cumprem papéis diferentes. A estimativa indica o tamanho e o sentido da associação; o erro-padrão indica a precisão dessa estimativa. Em termos simples, ele resume quanto a estimativa tenderia a variar se repetíssemos a coleta e o ajuste muitas vezes sob condições semelhantes. O valor de \(t\) é calculado por \(t=\widehat\beta/EP(\widehat\beta)\): para um mesmo coeficiente, maior erro-padrão implica menor precisão e menor evidência contra \(H_0:\beta=0\).

4.2 Modelo 1 — todas as medidas numéricas

O primeiro ajuste utiliza as três medidas disponíveis como preditoras e ainda não inclui a espécie:

\[ Y_i=\beta_0+\beta_1\,Sepal.Width_i+\beta_2\,Petal.Length_i+ \beta_3\,Petal.Width_i+\varepsilon_i. \]

modelo1 <- lm(
  Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width,
  data = dados
)

summary(modelo1)

## 
## Call:
## lm(formula = Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width, 
##     data = dados)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8282 -0.2199  0.0187  0.1971  0.8457 
## 
## Coefficients:
##              Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)    1.8560     0.2508    7.40      0.0000000000099 ***
## Sepal.Width    0.6508     0.0666    9.77 < 0.0000000000000002 ***
## Petal.Length   0.7091     0.0567   12.50 < 0.0000000000000002 ***
## Petal.Width   -0.5565     0.1275   -4.36      0.0000241287569 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.315 on 146 degrees of freedom
## Multiple R-squared:  0.859,  Adjusted R-squared:  0.856 
## F-statistic:  296 on 3 and 146 DF,  p-value: <0.0000000000000002

tabela_coeficientes(modelo1) |>
  knitr::kable(digits = 3) |>
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

termo	estimativa	erro-padrão	estatistica_t	p-valor	IC 95% inferior	IC 95% superior
(Intercept)	1.856	0.251	7.401	0	1.360	2.352
Sepal.Width	0.651	0.067	9.765	0	0.519	0.783
Petal.Length	0.709	0.057	12.502	0	0.597	0.821
Petal.Width	-0.556	0.128	-4.363	0	-0.809	-0.304

A equação ajustada é:

\[ \widehat{Sepal.Length}= 1,856 + 0,651\,Sepal.Width + 0,709\,Petal.Length - 0,556\,Petal.Width. \]

4.2.1 Comentários e observações sobre o Modelo 1

O modelo é globalmente significativo e apresenta \(R^2\) ajustado igual a 0,856, ou aproximadamente 85,6% da variabilidade explicada após a correção pelo número de parâmetros. Todas as variáveis têm p-valor inferior a 0,05.

Vamos isolar os erros-padrão dos três preditores para tornar a comparação explícita:

comparacao_ep_modelo1 <- tidy(modelo1) |>
  filter(term != "(Intercept)") |>
  transmute(
    variavel = term,
    estimativa = estimate,
    `erro-padrão` = std.error,
    `EP / |estimativa| (%)` = 100 * std.error / abs(estimate),
    `EP em relação ao menor EP` = std.error / min(std.error)
  )

comparacao_ep_modelo1 |>
  knitr::kable(digits = 3) |>
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

variavel	estimativa	erro-padrão	EP / |estimativa| (%)	EP em relação ao menor EP
Sepal.Width	0.651	0.067	10.240	1.175
Petal.Length	0.709	0.057	7.998	1.000
Petal.Width	-0.556	0.128	22.920	2.249

Dois resultados da própria saída do Modelo 1 merecem atenção:

1. o erro-padrão de Petal.Width é 0,128, contra 0,067 para Sepal.Width e 0,057 para Petal.Length. Portanto, ele é aproximadamente 1,91 vez o erro-padrão de Sepal.Width e 2,25 vezes o de Petal.Length;
2. Petal.Length apresenta coeficiente positivo (0,709), enquanto Petal.Width apresenta coeficiente negativo (-0,556), apesar da forte correlação positiva entre ambas.

Como todos os três preditores estão medidos em centímetros, a comparação direta dos erros-padrão é informativa neste exemplo. Mesmo assim, um erro-padrão maior não deve ser interpretado isoladamente: ele também depende da dispersão da variável, do tamanho da amostra e da relação desse preditor com os demais. A coluna percentual da tabela fornece uma segunda referência, comparando o erro-padrão com o tamanho absoluto da respectiva estimativa.

Essa inversão de sinal não é uma contradição matemática: na regressão múltipla, cada coeficiente representa uma associação parcial, mantendo os demais preditores constantes. Contudo, combinada com a correlação de aproximadamente 0,963, ela sugere que as duas medidas da pétala competem para explicar a mesma parcela da resposta. Para quantificar esse problema, calculamos o VIF.

4.2.2 O que é VIF?

VIF significa Variance Inflation Factor, ou fator de inflação da variância. Para o preditor \(X_j\),

\[ VIF_j=\frac{1}{1-R_j^2}, \]

em que \(R_j^2\) vem de uma regressão auxiliar de \(X_j\) contra todos os outros preditores. Se os demais preditores quase conseguem reconstruir \(X_j\), esse \(R_j^2\) fica próximo de 1 e o VIF cresce.

O VIF informa por quanto a variância estimada de \(\widehat\beta_j\) foi inflada pela relação linear entre os preditores. A raiz quadrada do VIF informa aproximadamente por quanto o erro-padrão foi inflado em comparação com uma situação ideal de preditores não correlacionados. A tolerância é \(1/VIF\); valores pequenos significam que resta pouca informação exclusiva naquele preditor.

vif_modelo1 <- performance::check_collinearity(modelo1) |>
  as.data.frame() |>
  transmute(
    variavel = Term,
    VIF = VIF,
    `raiz do VIF` = sqrt(VIF),
    tolerancia = Tolerance,
    R2_regressao_auxiliar = 1 - 1 / VIF
  )

vif_modelo1 |>
  rename(`R² da regressão auxiliar` = R2_regressao_auxiliar) |>
  knitr::kable(digits = 3) |>
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

variavel	VIF	raiz do VIF	tolerancia	R² da regressão auxiliar
Sepal.Width	1.271	1.127	0.787	0.213
Petal.Length	15.098	3.886	0.066	0.934
Petal.Width	14.234	3.773	0.070	0.930

Para Sepal.Width, o VIF é baixo (1,27). Para Petal.Length e Petal.Width, os VIF são 15,10 e 14,23. Isso significa que aproximadamente 93,4% da variação de Petal.Length e 93,0% da variação de Petal.Width podem ser explicadas pelos outros preditores do modelo. As raízes dos VIF, próximas de 3,89 e 3,77, evidenciam uma inflação importante da incerteza dos coeficientes.

Regras como “VIF maior que 5 exige atenção” e “VIF maior que 10 indica colinearidade forte” são referências práticas, não leis universais. Aqui, valores acima de 14, correlação de 0,963 e mudança de sinal contam uma história coerente: os efeitos individuais das duas medidas da pétala são difíceis de separar com estabilidade.

Leitura correta da multicolinearidade. O modelo ainda pode ajustar bem os dados. O principal prejuízo está na separação e interpretação dos efeitos individuais, nos erros-padrão e na estabilidade dos coeficientes.

4.3 Modelo 2 — retirada de Petal.Width

Como Petal.Length apresenta associação mais forte com a resposta e as duas medidas da pétala são muito correlacionadas, o segundo modelo retira Petal.Width:

\[ Y_i=\beta_0+\beta_1\,Sepal.Width_i+\beta_2\,Petal.Length_i+\varepsilon_i. \]

modelo2 <- lm(
  Sepal.Length ~ Sepal.Width + Petal.Length,
  data = dados
)

summary(modelo2)

## 
## Call:
## lm(formula = Sepal.Length ~ Sepal.Width + Petal.Length, data = dados)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9616 -0.2349  0.0008  0.2145  0.7856 
## 
## Coefficients:
##              Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)    2.2491     0.2480    9.07   0.0000000000000007 ***
## Sepal.Width    0.5955     0.0693    8.59   0.0000000000000116 ***
## Petal.Length   0.4719     0.0171   27.57 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.333 on 147 degrees of freedom
## Multiple R-squared:  0.84,   Adjusted R-squared:  0.838 
## F-statistic:  386 on 2 and 147 DF,  p-value: <0.0000000000000002

tabela_coeficientes(modelo2) |>
  knitr::kable(digits = 3) |>
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

termo	estimativa	erro-padrão	estatistica_t	p-valor	IC 95% inferior	IC 95% superior
(Intercept)	2.249	0.248	9.07	0	1.759	2.739
Sepal.Width	0.596	0.069	8.59	0	0.459	0.733
Petal.Length	0.472	0.017	27.57	0	0.438	0.506

A equação ajustada é:

\[ \widehat{Sepal.Length}= 2,249 + 0,596\,Sepal.Width + 0,472\,Petal.Length. \]

4.3.1 Comentários e interpretação do Modelo 2

O modelo permanece globalmente significativo e explica aproximadamente 83,8% da variabilidade do comprimento da sépala, segundo o \(R^2\) ajustado.

• Mantendo o comprimento da pétala constante, o aumento de 1 cm em Sepal.Width está associado a um aumento médio estimado de 0,596 cm em Sepal.Length.
• Mantendo a largura da sépala constante, o aumento de 1 cm em Petal.Length está associado a um aumento médio estimado de 0,472 cm em Sepal.Length.

As duas associações são positivas, estatisticamente significativas e simples de comunicar. A retirada de Petal.Width reduz a colinearidade e custa cerca de 1,8 ponto percentual de \(R^2\) ajustado em relação ao Modelo 1.

4.4 Modelo 3 — retirada de Petal.Length

Para avaliar a alternativa oposta, retiramos Petal.Length e mantemos Petal.Width:

\[ Y_i=\beta_0+\beta_1\,Sepal.Width_i+\beta_2\,Petal.Width_i+\varepsilon_i. \]

modelo3 <- lm(
  Sepal.Length ~ Sepal.Width + Petal.Width,
  data = dados
)

summary(modelo3)

## 
## Call:
## lm(formula = Sepal.Length ~ Sepal.Width + Petal.Width, data = dados)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.208 -0.229 -0.045  0.227  1.181 
## 
## Coefficients:
##             Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)   3.4573     0.3092   11.18 < 0.0000000000000002 ***
## Sepal.Width   0.3991     0.0911    4.38             0.000022 ***
## Petal.Width   0.9721     0.0521   18.66 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.451 on 147 degrees of freedom
## Multiple R-squared:  0.707,  Adjusted R-squared:  0.703 
## F-statistic:  178 on 2 and 147 DF,  p-value: <0.0000000000000002

tabela_coeficientes(modelo3) |>
  knitr::kable(digits = 3) |>
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

termo	estimativa	erro-padrão	estatistica_t	p-valor	IC 95% inferior	IC 95% superior
(Intercept)	3.457	0.309	11.18	0	2.846	4.068
Sepal.Width	0.399	0.091	4.38	0	0.219	0.579
Petal.Width	0.972	0.052	18.66	0	0.869	1.075

A equação ajustada é:

\[ \widehat{Sepal.Length}= 3,457 + 0,399\,Sepal.Width + 0,972\,Petal.Width. \]

4.4.1 Comentários e interpretação do Modelo 3

As duas variáveis também são significativas, mas o \(R^2\) ajustado cai para 0,703. Mantendo Petal.Width constante, o aumento de 1 cm em Sepal.Width está associado a 0,399 cm adicionais no comprimento médio da sépala. Mantendo Sepal.Width constante, 1 cm adicional em Petal.Width está associado a 0,972 cm adicionais na resposta.

O ajuste inferior ao do Modelo 2 é coerente com a exploração: Petal.Width tem correlação total um pouco menor com Sepal.Length que Petal.Length. Entre os dois modelos que eliminam uma das medidas colineares, o Modelo 2 preserva mais informação.

4.5 Preparação para os Modelos 4 e 5 — variável Species

A espécie possui três categorias. O R usa setosa como referência e cria duas variáveis indicadoras:

Espécie	\(Z_1\): versicolor	\(Z_2\): virginica
setosa	0	0
versicolor	1	0
virginica	0	1

levels(dados$Species)

## [1] "setosa"     "versicolor" "virginica"

contrasts(dados$Species)

##            versicolor virginica
## setosa              0         0
## versicolor          1         0
## virginica           0         1

amostra_dummies <- dados[c(1, 51, 101), ] |>
  bind_cols(
    as.data.frame(model.matrix(~ Species, data = dados[c(1, 51, 101), ])[, -1])
  )

amostra_dummies |>
  select(id, Species, Speciesversicolor, Speciesvirginica) |>
  knitr::kable() |>
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

	id	Species	Speciesversicolor	Speciesvirginica
1	1	setosa	0	0
51	51	versicolor	1	0
101	101	virginica	0	1

Assim, dois coeficientes bastam para representar três categorias sem criar dependência linear perfeita com o intercepto.

4.6 Modelo 4 — medidas selecionadas e espécie

O quarto modelo acrescenta Species ao Modelo 2:

\[ Y_i=\beta_0+\beta_1\,Sepal.Width_i+\beta_2\,Petal.Length_i+ \beta_3 Z_{1i}+\beta_4 Z_{2i}+\varepsilon_i. \]

modelo4 <- lm(
  Sepal.Length ~ Sepal.Width + Petal.Length + Species,
  data = dados
)

summary(modelo4)

## 
## Call:
## lm(formula = Sepal.Length ~ Sepal.Width + Petal.Length + Species, 
##     data = dados)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8216 -0.2053  0.0064  0.2265  0.7500 
## 
## Coefficients:
##                   Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)         2.3904     0.2623    9.11  0.00000000000000059 ***
## Sepal.Width         0.4322     0.0814    5.31  0.00000040259819018 ***
## Petal.Length        0.7756     0.0642   12.07 < 0.0000000000000002 ***
## Speciesversicolor  -0.9558     0.2152   -4.44  0.00001759998870080 ***
## Speciesvirginica   -1.3941     0.2857   -4.88  0.00000275961767052 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.31 on 145 degrees of freedom
## Multiple R-squared:  0.863,  Adjusted R-squared:  0.86 
## F-statistic:  229 on 4 and 145 DF,  p-value: <0.0000000000000002

tabela_coeficientes(modelo4) |>
  knitr::kable(digits = 3) |>
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

termo	estimativa	erro-padrão	estatistica_t	p-valor	IC 95% inferior	IC 95% superior
(Intercept)	2.390	0.262	9.114	0	1.872	2.909
Sepal.Width	0.432	0.081	5.310	0	0.271	0.593
Petal.Length	0.776	0.064	12.073	0	0.649	0.903
Speciesversicolor	-0.956	0.215	-4.442	0	-1.381	-0.530
Speciesvirginica	-1.394	0.286	-4.880	0	-1.959	-0.830

Com setosa como referência, a equação ajustada é:

\[ \widehat{Sepal.Length}= 2,390 + 0,432\,Sepal.Width + 0,776\,Petal.Length - 0,956\,Z_1 - 1,394\,Z_2. \]

4.6.1 Comentários e observações sobre o Modelo 4

O Modelo 4 apresenta o maior \(R^2\) ajustado entre os cinco ajustes (0,860) e todos os termos são significativos. Entretanto, os coeficientes de versicolor e virginica ficam negativos em relação a setosa, o que parece estranho diante dos gráficos, nos quais essas espécies têm sépalas maiores.

Para compreender esse resultado, observe que Species e Petal.Length estão fortemente associados. Os coeficientes de espécie no Modelo 4 não comparam as médias brutas das espécies; eles comparam flores de espécies diferentes que possuam os mesmos valores de Sepal.Width e Petal.Length. Como há pouca sobreposição de Petal.Length entre setosa e as outras espécies, essa interpretação condicional é frágil e próxima de uma extrapolação.

vif_modelo4 <- performance::check_collinearity(modelo4) |>
  as.data.frame() |>
  transmute(
    variavel = Term,
    `VIF/GVIF bruto` = VIF,
    `graus de liberdade` = ifelse(Term == "Species", nlevels(dados$Species) - 1, 1),
    `fator comparável de inflação do EP` = VIF^(1 / (2 * `graus de liberdade`)),
    tolerancia = Tolerance
  )

vif_modelo4 |>
  knitr::kable(digits = 3) |>
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

variavel	VIF/GVIF bruto	graus de liberdade	fator comparável de inflação do EP	tolerancia
Sepal.Width	1.947	1	1.395	0.514
Petal.Length	19.899	1	4.461	0.050
Species	27.112	2	2.282	0.037

Para Petal.Length, o VIF é 19,90 e o fator de inflação do erro-padrão é 4,46. Para Species, que possui dois graus de liberdade, o valor bruto é um GVIF de 27,11; para compará-lo com variáveis de um grau de liberdade, usamos o fator ajustado GVIF^(1/(2×gl)), igual a 2,28. Portanto, o Modelo 4 oferece melhor ajuste interno, mas seus coeficientes condicionais — sobretudo os contrastes de espécie, apoiados por pouca sobreposição de Petal.Length — exigem cautela e não devem ser tratados como efeitos estáveis ou causais.

4.7 Modelo 5 — largura da sépala e espécie

Para obter contrastes de espécie mais diretamente compatíveis com a exploração gráfica, retiramos Petal.Length:

\[ Y_i=\beta_0+\beta_1\,Sepal.Width_i+\beta_2 Z_{1i}+\beta_3 Z_{2i}+\varepsilon_i. \]

modelo5 <- lm(
  Sepal.Length ~ Sepal.Width + Species,
  data = dados
)

summary(modelo5)

## 
## Call:
## lm(formula = Sepal.Length ~ Sepal.Width + Species, data = dados)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.3071 -0.2571 -0.0533  0.1954  1.4125 
## 
## Coefficients:
##                   Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)          2.251      0.370    6.09      0.0000000095681 ***
## Sepal.Width          0.804      0.106    7.56      0.0000000000042 ***
## Speciesversicolor    1.459      0.112   13.01 < 0.0000000000000002 ***
## Speciesvirginica     1.947      0.100   19.47 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.438 on 146 degrees of freedom
## Multiple R-squared:  0.726,  Adjusted R-squared:  0.72 
## F-statistic:  129 on 3 and 146 DF,  p-value: <0.0000000000000002

tabela_coeficientes(modelo5) |>
  knitr::kable(digits = 3) |>
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

termo	estimativa	erro-padrão	estatistica_t	p-valor	IC 95% inferior	IC 95% superior
(Intercept)	2.251	0.370	6.089	0	1.521	2.982
Sepal.Width	0.804	0.106	7.557	0	0.593	1.014
Speciesversicolor	1.459	0.112	13.012	0	1.237	1.680
Speciesvirginica	1.947	0.100	19.465	0	1.749	2.144

Com setosa como referência, a equação ajustada é:

\[ \widehat{Sepal.Length}= 2,251 + 0,804\,Sepal.Width + 1,459\,Z_1 + 1,947\,Z_2. \]

4.7.1 Comentários e interpretação do Modelo 5

Agora os coeficientes de espécie acompanham as diferenças observadas nos gráficos:

• controlando Sepal.Width, flores versicolor apresentam comprimento médio de sépala 1,459 cm maior que flores setosa;
• controlando Sepal.Width, flores virginica apresentam comprimento médio de sépala 1,947 cm maior que flores setosa.

As interpretações são mais intuitivas, mas o \(R^2\) ajustado cai para 0,720, bem abaixo do Modelo 2. Portanto, retirar Petal.Length melhora a leitura dos contrastes de espécie, porém elimina uma variável quantitativa muito informativa.

4.8 Alterando a categoria de referência

Mudar a referência altera a forma dos coeficientes, mas não altera resíduos, valores ajustados, \(R^2\), AIC ou BIC.

dados_ref_virginica <- dados |>
  mutate(Species = relevel(Species, ref = "virginica"))

modelo5_ref_virginica <- lm(
  Sepal.Length ~ Sepal.Width + Species,
  data = dados_ref_virginica
)

tabela_coeficientes(modelo5_ref_virginica) |>
  knitr::kable(digits = 3) |>
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

termo	estimativa	erro-padrão	estatistica_t	p-valor	IC 95% inferior	IC 95% superior
(Intercept)	4.198	0.322	13.027	0	3.561	4.835
Sepal.Width	0.804	0.106	7.557	0	0.593	1.014
Speciessetosa	-1.947	0.100	-19.465	0	-2.144	-1.749
Speciesversicolor	-0.488	0.090	-5.409	0	-0.666	-0.310

Com virginica como referência, os coeficientes negativos de setosa e versicolor expressam quanto essas espécies têm sépalas menores, em média condicional, do que virginica. O exemplo mostra que mudar a referência muda as comparações expressas pelos coeficientes, mas não muda o ajuste do modelo.

4.9 Seleção do modelo

Antes de examinar resíduos ou fazer previsões, comparamos os modelos e identificamos um candidato de melhor ajuste. A tabela abaixo reúne os cinco ajustes. Valores maiores de \(R^2\) ajustado e menores de RQME residual, AIC e BIC são preferíveis.

• \(R^2\) ajustado: parcela da variação explicada, com penalização pelo número de parâmetros;
• RQME residual: tamanho típico do erro de ajuste, em centímetros;
• AIC e BIC: critérios que equilibram ajuste e complexidade; valores menores são preferíveis, e o BIC penaliza mais modelos complexos.

comparacao_ajuste <- imap_dfr(modelos, function(ajuste, nome) {
  g <- glance(ajuste)
  tibble(
    modelo = nome,
    parametros = length(coef(ajuste)),
    R2 = g$r.squared,
    R2_ajustado = g$adj.r.squared,
    RQME_residual = sqrt(mean(residuals(ajuste)^2)),
    AIC = AIC(ajuste),
    BIC = BIC(ajuste)
  )
})

knitr::kable(comparacao_ajuste, digits = 3) |>
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

modelo	parametros	R2	R2_ajustado	RQME_residual	AIC	BIC
Modelo 1 — todas as medidas	4	0.859	0.856	0.310	84.64	99.70
Modelo 2 — sem Petal.Width	3	0.840	0.838	0.330	101.03	113.07
Modelo 3 — sem Petal.Length	3	0.707	0.703	0.447	191.82	203.86
Modelo 4 — medidas + espécie	5	0.863	0.860	0.305	81.58	99.64
Modelo 5 — largura + espécie	4	0.726	0.720	0.432	183.94	198.99

O Modelo 4 lidera os critérios de ajuste interno entre os cinco modelos. O Modelo 1 fica próximo em ajuste, mas mantém a multicolinearidade entre as duas medidas da pétala. O Modelo 2 é mais simples, mas perde ajuste; os Modelos 3 e 5 perdem informação importante ao excluir Petal.Length.

Como complemento, os testes F abaixo verificam se um termo específico acrescenta informação quando comparamos modelos aninhados. Cada comparação é apresentada separadamente para que fique claro qual variável está sendo avaliada.

4.9.1 Comparação 1 — Modelo 2 versus Modelo 1: contribuição de Petal.Width

O Modelo 1 contém Sepal.Width, Petal.Length e Petal.Width; o Modelo 2 retira apenas Petal.Width. Este teste verifica se essa variável acrescenta ajuste depois de considerar as outras duas medidas.

anova(modelo2, modelo1)

4.9.2 Comparação 2 — Modelo 3 versus Modelo 1: contribuição de Petal.Length

O Modelo 3 retira Petal.Length do Modelo 1. Assim, o teste avalia a contribuição dessa medida quando Sepal.Width e Petal.Width já estão no modelo.

anova(modelo3, modelo1)

4.9.3 Comparação 3 — Modelo 2 versus Modelo 4: contribuição conjunta de Species

O Modelo 4 acrescenta Species ao Modelo 2. Este é o teste que verifica diretamente se representar as diferenças entre espécies melhora o ajuste.

teste_m2_m4 <- anova(modelo2, modelo4)
teste_m2_m4

4.9.4 Comparação 4 — Modelo 5 versus Modelo 4: contribuição de Petal.Length na presença de Species

O Modelo 5 retira Petal.Length do Modelo 4. O teste mostra se essa medida continua trazendo informação mesmo depois de controlar a espécie.

anova(modelo5, modelo4)

A tabela identifica o candidato de melhor ajuste interno; os testes F verificam se ganhos específicos são compatíveis com os dados. Em especial, o teste anova(modelo2, modelo4) avalia diretamente a passagem do Modelo 2 para o Modelo 4: a inclusão conjunta de Species melhora significativamente o ajuste (p-valor 0,000012). O teste F responde a uma pergunta pontual — “este termo melhora este modelo reduzido?” — e só pode ser usado para modelos aninhados. AIC, BIC e \(R^2\) ajustado comparam o equilíbrio geral entre ajuste e complexidade.

Candidato de melhor ajuste interno: Modelo 4 — Sepal.Length ~ Sepal.Width + Petal.Length + Species. Ele apresenta o maior \(R^2\) ajustado (0,860), o menor RQME residual (0,305), o menor AIC (81,57) e o menor BIC (99,64) entre os cinco ajustes. Como a colinearidade é forte, a decisão final depende também do diagnóstico e do objetivo: ajuste interno ou interpretação mais estável. Esta comparação não demonstra desempenho em novas flores; isso exigiria uma base de teste.

5 Diagnóstico completo dos resíduos

O Modelo 4, candidato de melhor ajuste na tabela anterior, será usado para a análise detalhada dos pressupostos. O diagnóstico verificará se sua forma linear é compatível com os dados antes da decisão final.

modelo_diagnostico <- modelo4

diag <- augment(modelo_diagnostico, data = dados) |>
  mutate(
    indice = row_number(),
    cook_limite = 4 / nrow(dados),
    leverage_limite = 2 * length(coef(modelo_diagnostico)) / nrow(dados)
  )

5.1 O que verificaremos antes dos gráficos

O diagnóstico não procura um “gráfico perfeito”; ele verifica se há indícios fortes de que o modelo linear escolhido é inadequado. Antes de olhar a figura, estabeleça as perguntas de cada painel:

• Linearidade: os resíduos oscilam sem curvatura em torno de zero?
• Homocedasticidade: a dispersão dos resíduos se mantém semelhante ao longo dos valores ajustados?
• Normalidade aproximada: a distribuição dos resíduos é compatível com uma normal, especialmente nas caudas?
• Alavancagem e influência: existe alguma flor com combinação incomum de preditores e capacidade de alterar substancialmente o ajuste?

Essas verificações se complementam. Um único ponto fora de um limite não invalida o modelo; o interesse está em padrões persistentes e em sua relevância para a pergunta da análise.

5.2 Quatro gráficos fundamentais

p1 <- ggplot(diag, aes(.fitted, .resid)) +
  geom_hline(yintercept = 0, color = cinza, linetype = 2) +
  geom_point(color = azul, alpha = 0.7) +
  geom_smooth(se = FALSE, method = "loess", color = rosa, linewidth = 0.85) +
  labs(title = "Resíduos × ajustados", subtitle = "Procure curvatura e formato de funil",
       x = "Valor ajustado", y = "Resíduo")

set.seed(20260715)
n_diag <- nrow(diag)
quantis_teoricos <- qnorm(ppoints(n_diag))
simulacoes_normais <- replicate(2000, sort(rnorm(n_diag)))
envelope_qq <- tibble(
  quantil_teorico = quantis_teoricos,
  residuo_observado = sort(diag$.std.resid),
  limite_inferior = apply(simulacoes_normais, 1, quantile, probs = 0.025),
  limite_superior = apply(simulacoes_normais, 1, quantile, probs = 0.975)
)

p2 <- ggplot(envelope_qq, aes(quantil_teorico, residuo_observado)) +
  geom_ribbon(
    aes(ymin = limite_inferior, ymax = limite_superior),
    fill = dourado, alpha = 0.22
  ) +
  geom_abline(intercept = 0, slope = 1, color = rosa, linewidth = 0.85) +
  geom_point(color = azul, alpha = 0.72) +
  labs(title = "QQ-plot com envelope de 95%", subtitle = "Pontos fora da faixa sugerem desvio de normalidade",
       x = "Quantis teóricos", y = "Resíduo padronizado")

p3 <- ggplot(diag, aes(.fitted, sqrt(abs(.std.resid)))) +
  geom_point(color = azul, alpha = 0.7) +
  geom_smooth(se = FALSE, method = "loess", color = rosa, linewidth = 0.85) +
  labs(title = "Escala–localização", subtitle = "Uma faixa aproximadamente horizontal é desejável",
       x = "Valor ajustado", y = expression(sqrt("|resíduo padronizado|")))

p4 <- ggplot(diag, aes(.hat, .std.resid, size = .cooksd)) +
  geom_hline(yintercept = c(-2, 2), color = cinza, linetype = 2) +
  geom_vline(xintercept = unique(diag$leverage_limite), color = dourado, linetype = 2) +
  geom_point(color = azul, alpha = 0.62) +
  scale_size_continuous(range = c(1.2, 6), guide = "none") +
  labs(title = "Alavancagem e influência", subtitle = "Linha laranja: alavancagem alta; tamanho: distância de Cook",
       x = "Alavancagem", y = "Resíduo padronizado")

(p1 + p2) / (p3 + p4)

Diagnóstico do Modelo 4 com as 150 observações. Linhas de referência servem como guias, não como decisões mecânicas.

5.2.1 Linearidade

O primeiro gráfico coloca no eixo horizontal os valores previstos pelo modelo e, no eixo vertical, os resíduos \(e_i=y_i-\widehat y_i\). Um resíduo positivo indica que o modelo subestimou o comprimento da sépala; um resíduo negativo indica que o modelo superestimou.

O gráfico responde principalmente à pergunta: a forma linear escolhida consegue representar a média da resposta? O comportamento desejável é uma nuvem distribuída aleatoriamente ao redor de zero, com a linha suavizada aproximadamente horizontal. Uma forma em U ou U invertido sugeriria curvatura não representada; ondas poderiam indicar termos omitidos; grupos separados poderiam revelar que alguma estrutura categórica ainda não foi modelada.

No Modelo 4, a linha suavizada permanece próxima de zero e não aparece uma curvatura persistente. Há resíduos positivos e negativos ao longo de toda a faixa de valores ajustados. Portanto, o gráfico não mostra evidência visual forte de que a média precise de um termo não linear. Isso é compatível com linearidade, mas não constitui uma prova definitiva.

5.2.2 Homocedasticidade

O terceiro gráfico, escala–localização, retira o sinal dos resíduos padronizados e destaca seu tamanho. Ele responde à pergunta: a dispersão dos erros permanece aproximadamente constante para valores previstos pequenos e grandes? Uma faixa de largura semelhante e uma linha suavizada horizontal são desejáveis. Um funil que se abre indicaria variância crescente; um funil que se fecha indicaria variância decrescente.

No Modelo 4, não se observa um funil forte, e a linha suavizada varia pouco ao longo do eixo horizontal. O gráfico de resíduos × ajustados conduz à mesma conclusão visual: a amplitude vertical não aumenta sistematicamente com o valor previsto. Assim, a inspeção gráfica não sugere heterocedasticidade importante.

5.2.3 Normalidade aproximada

O QQ-plot compara os quantis dos resíduos padronizados com os quantis esperados de uma distribuição normal. A faixa dourada é um envelope simulado heurístico de 95%: ela mostra a variação esperada de uma amostra normal independente do mesmo tamanho. Resíduos padronizados de uma regressão não são observações independentes perfeitas; por isso, o envelope é apenas uma referência visual. Pontos sistematicamente fora da faixa, sobretudo nas caudas, sugerem desvio de normalidade; poucos pontos isolados podem ocorrer ao acaso.

No Modelo 4, os pontos acompanham bem a reta, inclusive nas caudas, sem uma forma sistemática de S. O QQ-plot avalia a distribuição dos resíduos, e não exige que Sepal.Length ou os preditores sejam normalmente distribuídos. Essa normalidade é especialmente relevante para a qualidade dos testes, p-valores e intervalos de confiança em amostras pequenas. O envelope não avalia alavancagem ou influência: uma observação pode estar dentro do envelope de normalidade e, ainda assim, ter alavancagem alta.

teste_shapiro <- shapiro.test(residuals(modelo_diagnostico))
tibble(
  teste = "Shapiro–Wilk",
  W = unname(teste_shapiro$statistic),
  `p-valor` = teste_shapiro$p.value
) |>
  mutate(across(where(is.numeric), ~ round(.x, 4)))

O teste de Shapiro–Wilk apresentou p-valor 0,9048. Como esse valor é maior que 0,05, não rejeitamos a hipótese de normalidade. O teste e o QQ-plot apontam na mesma direção. Ainda assim, com amostras grandes, testes podem detectar desvios irrelevantes; com amostras pequenas, podem ter baixa potência. Por isso, o gráfico e o teste devem ser lidos em conjunto.

5.2.4 Alavancagem e influência

O quarto gráfico combina três informações. A alavancagem, no eixo horizontal, mede quão incomum é a combinação de preditores de uma flor. A linha laranja é apenas um limite de triagem para alavancagem alta; ultrapassá-la significa que a flor está em uma região pouco comum dos preditores, e não que ela seja uma observação errada ou que deva ser removida. O resíduo padronizado, no eixo vertical, mede quão distante sua resposta ficou do valor previsto. O tamanho do ponto representa a distância de Cook, que resume quanto o ajuste completo pode mudar quando aquela observação é retirada.

Uma observação com alavancagem alta não é automaticamente problemática: ela pode seguir perfeitamente a tendência. Uma observação com resíduo grande também pode ter pouca influência se estiver em uma região muito povoada dos preditores. A maior preocupação ocorre quando alavancagem alta e resíduo grande aparecem juntos, produzindo distância de Cook elevada.

Usando apenas referências de triagem, 7 observações ultrapassam $2p/n=0,067 e 7 ultrapassam $4/n=0,027 para a distância de Cook. A maior distância de Cook é 0,067, muito abaixo da referência ampla de 1. Isso sugere pontos que merecem inspeção, mas não uma observação isolada dominando completamente o Modelo 4.

5.2.5 Independência

A independência não pode ser confirmada por um gráfico se o desenho da coleta for desconhecido. Se várias flores vieram da mesma planta, canteiro, região ou lote, talvez seja necessário um modelo hierárquico. Se as linhas estiverem ordenadas no tempo, avalie autocorrelação.

ggplot(diag, aes(indice, .resid)) +
  geom_hline(yintercept = 0, color = cinza, linetype = 2) +
  geom_line(color = "#9FB4D1", linewidth = 0.55) +
  geom_point(color = azul, size = 1.5) +
  labs(
    title = "Resíduos na ordem das observações",
    subtitle = "Só é informativo sobre dependência se a ordem tiver significado no desenho",
    x = "Ordem", y = "Resíduo"
  )

5.2.6 Investigando as observações sinalizadas

influentes <- augment(modelo_diagnostico, data = dados) |>
  mutate(
    id = dados$id,
    limiar_cook = 4 / nrow(dados),
    sinalizada = .cooksd > limiar_cook | abs(.std.resid) > 2
  ) |>
  arrange(desc(.cooksd)) |>
  transmute(
    id,
    especie = Species,
    comprimento_sepala = Sepal.Length,
    largura_sepala = Sepal.Width,
    comprimento_petala = Petal.Length,
    valor_ajustado = .fitted,
    residuo = .resid,
    residuo_padronizado = .std.resid,
    alavancagem = .hat,
    distancia_cook = .cooksd,
    sinalizada
  )

DT::datatable(
  head(influentes, 12),
  rownames = FALSE,
  options = list(pageLength = 6, scrollX = TRUE),
  caption = "Observações com maior distância de Cook",
  colnames = c(
    "ID", "Espécie", "Comprimento da sépala", "Largura da sépala",
    "Comprimento da pétala", "Valor ajustado", "Resíduo",
    "Resíduo padronizado", "Alavancagem", "Distância de Cook", "Sinalizada"
  )
) |>
  formatRound(columns = c("comprimento_sepala", "largura_sepala", "comprimento_petala",
                          "valor_ajustado", "residuo", "residuo_padronizado",
                          "alavancagem", "distancia_cook"), digits = 3)

Distância de Cook acima de \(4/n\) é um sinal para investigar, não uma ordem para excluir. Primeiro verifique erro de digitação, compatibilidade com a população e sensibilidade das conclusões com e sem o ponto.

Análise de sensibilidade, não seleção de modelo. Os pontos abaixo são retirados apenas para verificar a estabilidade das conclusões. A comparação não deve ser usada para escolher um modelo com \(R^2\) maior nem para justificar exclusões automáticas.

ids_sensiveis <- influentes |>
  filter(sinalizada) |>
  pull(id)

modelo_sem_sensiveis <- lm(
  formula(modelo_diagnostico),
  data = dados |> filter(!id %in% ids_sensiveis)
)

metricas_sensibilidade <- tibble(
  metrica = c("R²", "R² ajustado", "Erro-padrão residual"),
  `Todos os dados` = c(
    summary(modelo_diagnostico)$r.squared,
    summary(modelo_diagnostico)$adj.r.squared,
    summary(modelo_diagnostico)$sigma
  ),
  `Sem pontos sinalizados` = c(
    summary(modelo_sem_sensiveis)$r.squared,
    summary(modelo_sem_sensiveis)$adj.r.squared,
    summary(modelo_sem_sensiveis)$sigma
  )
)

metricas_sensibilidade |>
  knitr::kable(digits = 3, caption = "Métricas de ajuste na análise de sensibilidade") |>
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Métricas de ajuste na análise de sensibilidade

metrica	Todos os dados	Sem pontos sinalizados
R²	0.863	0.898
R² ajustado	0.860	0.895
Erro-padrão residual	0.310	0.269

comparacao_sensibilidade <- bind_rows(
  tidy(modelo_diagnostico) |> mutate(ajuste = "Todos os dados"),
  tidy(modelo_sem_sensiveis) |> mutate(ajuste = "Sem pontos sinalizados")
) |>
  select(ajuste, term, estimate, std.error) |>
  pivot_wider(names_from = ajuste, values_from = c(estimate, std.error))

knitr::kable(comparacao_sensibilidade, digits = 3) |>
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

term	estimate_Todos os dados	estimate_Sem pontos sinalizados	std.error_Todos os dados	std.error_Sem pontos sinalizados
(Intercept)	2.390	2.358	0.262	0.233
Sepal.Width	0.432	0.423	0.081	0.072
Petal.Length	0.776	0.809	0.064	0.059
Speciesversicolor	-0.956	-0.981	0.215	0.193
Speciesvirginica	-1.394	-1.502	0.286	0.259

Foram sinalizadas 8 observações. Ao retirá-las apenas para esta análise de sensibilidade, o \(R^2\) ajustado passa de 0,860 para 0,895 e o erro residual passa de 0,310 para 0,269. Essa aparente “melhora” é esperada quando removemos pontos com resíduos maiores; ela não prova que o modelo sem essas flores é melhor ou que elas devam ser excluídas.

Os coeficientes permanecem próximos, como mostra a tabela. Portanto, o Modelo 4 com as 150 flores continua sendo o resultado principal. Uma observação só deveria ser removida se houver evidência independente de erro de registro, de que ela não pertence à população estudada ou uma regra definida antes da análise — nunca apenas por cruzar a linha laranja, ter resíduo grande ou ficar dentro/fora do envelope.

6 Previsão e incerteza

Com o Modelo 4 selecionado, podemos estimar o comprimento da sépala para novos valores de Sepal.Width, Petal.Length e Species. Os nomes e tipos das colunas no novo data.frame precisam coincidir com os usados no ajuste.

modelo_final <- modelo4

novas_flores <- tibble(
  cenario = c("Flor A", "Flor B", "Flor C"),
  Sepal.Width = c(2.5, 3.0, 3.5),
  Petal.Length = c(1.5, 4.0, 5.5),
  Species = factor(
    c("setosa", "versicolor", "virginica"),
    levels = levels(dados$Species)
  )
)

ic_media <- predict(modelo_final, newdata = novas_flores,
                    interval = "confidence", level = 0.95)
ip_individual <- predict(modelo_final, newdata = novas_flores,
                         interval = "prediction", level = 0.95)

bind_rows(
  bind_cols(novas_flores, as_tibble(ic_media)) |>
    mutate(intervalo = "IC 95% da média"),
  bind_cols(novas_flores, as_tibble(ip_individual)) |>
    mutate(intervalo = "IP 95% individual")
) |>
  select(cenario, Species, Sepal.Width, Petal.Length, intervalo,
         estimativa = fit, inferior = lwr, superior = upr) |>
  knitr::kable(digits = 3) |>
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

cenario	Species	Sepal.Width	Petal.Length	intervalo	estimativa	inferior	superior
Flor A	setosa	2.5	1.5	IC 95% da média	4.634	4.460	4.809
Flor B	versicolor	3.0	4.0	IC 95% da média	5.834	5.729	5.938
Flor C	virginica	3.5	5.5	IC 95% da média	6.775	6.652	6.898
Flor A	setosa	2.5	1.5	IP 95% individual	4.634	3.997	5.272
Flor B	versicolor	3.0	4.0	IP 95% individual	5.834	5.212	6.456
Flor C	virginica	3.5	5.5	IP 95% individual	6.775	6.149	7.401

• Intervalo de confiança da média: captura a incerteza sobre o comprimento médio esperado de todas as flores com aquele conjunto de preditores.
• Intervalo de previsão: captura a incerteza para uma nova flor individual; incorpora a variação residual e, por isso, é sempre mais amplo para os mesmos valores dos preditores.

Para uma previsão individual, use o intervalo de previsão. Usar o intervalo de confiança da média nesse caso produz uma faixa estreita demais e aumenta a chance de não conter o valor individual observado.

Dica para uma análise voltada à previsão. Em um projeto preditivo, separe os dados logo no início em duas bases: treinamento e teste. Use somente a base de treinamento para explorar os dados, ajustar os modelos, comparar especificações e diagnosticar resíduos — ou seja, para realizar todo o percurso apresentado até aqui. Guarde a base de teste sem usá-la nesse processo. Ao final, aplique predict() ao conjunto de teste e compare as previsões com os valores verdadeiros de Sepal.Length, por exemplo com RQME, MAE e gráficos de observado versus previsto. Essa avaliação mostra como o modelo se comporta em flores que não participaram do ajuste. Nesta aula, a base completa foi usada para privilegiar a explicação e a inferência; a divisão é uma extensão recomendada quando a finalidade principal for previsão.

7 Conclusões finais

7.1 Síntese da exploração

A base contém 150 flores, igualmente distribuídas entre setosa, versicolor e virginica, sem valores ausentes. As medidas das pétalas são fortemente correlacionadas entre si (\(r\approx0{,}963\)), o que explica a multicolinearidade observada quando ambas entram simultaneamente. A exploração também mostra separação marcante entre espécies, principalmente nas medidas da pétala; por isso, Species não pode ser ignorada sem uma justificativa substantiva.

7.2 Decisão e interpretação

Após os diagnósticos, o Modelo 4 permanece o melhor ajuste interno entre os cinco modelos: os gráficos não revelam um motivo claro para descartá-lo, e ele lidera os critérios de ajuste. Porém, ele não é automaticamente o melhor modelo para toda finalidade. Para descrição global interna, ele é a escolha principal; para ensinar e interpretar efeitos contínuos com menor colinearidade, o Modelo 2 é uma alternativa mais estável e mais simples. Não há base de teste nesta análise, portanto não é possível afirmar qual dos cinco é melhor para prever novas flores.

No Modelo 4:

• mantendo Petal.Length e Species constantes, 1 cm adicional em Sepal.Width está associado a um aumento médio de 0,432 cm em Sepal.Length;
• mantendo Sepal.Width e Species constantes, 1 cm adicional em Petal.Length está associado a um aumento médio de 0,776 cm em Sepal.Length;
• os coeficientes de versicolor e virginica são negativos em relação a setosa depois de fixar as duas medidas numéricas. Eles não contradizem os gráficos de médias brutas: respondem a uma comparação condicional diferente e apoiada por pouca sobreposição de comprimento de pétala entre as espécies.

tabela_coeficientes(modelo4) |>
  knitr::kable(digits = 3) |>
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

termo	estimativa	erro-padrão	estatistica_t	p-valor	IC 95% inferior	IC 95% superior
(Intercept)	2.390	0.262	9.114	0	1.872	2.909
Sepal.Width	0.432	0.081	5.310	0	0.271	0.593
Petal.Length	0.776	0.064	12.073	0	0.649	0.903
Speciesversicolor	-0.956	0.215	-4.442	0	-1.381	-0.530
Speciesvirginica	-1.394	0.286	-4.880	0	-1.959	-0.830

7.3 Verificação final dos pressupostos

Os gráficos de resíduos e escala–localização não sugerem mudança forte de variância. O teste de Shapiro–Wilk apresentou p-valor 0,9048, compatível com a normalidade aproximada indicada pelo QQ-plot com envelope. Os gráficos de alavancagem e distância de Cook foram examinados em conjunto e não revelaram um motivo isolado para descartar o ajuste.

A independência não pode ser estabelecida apenas com a base. Se houver várias flores provenientes da mesma planta, canteiro ou local, seria necessário representar essa estrutura no modelo. As observações sinalizadas por influência também devem ser investigadas, não eliminadas automaticamente.

7.4 Conclusão substantiva

O comprimento da sépala está positivamente associado à largura da sépala e ao comprimento da pétala, depois de controlar a espécie. A espécie organiza grande parte da estrutura morfológica da base e altera a interpretação das relações entre as medidas. Entre os cinco modelos estudados, o Modelo 4 oferece o melhor ajuste interno; o Modelo 2 oferece a leitura mais simples dos efeitos contínuos. Os contrastes de espécie do Modelo 4 precisam ser comunicados como associações condicionais, sob forte colinearidade.

As conclusões são associativas, não causais. A base Iris é observacional e não contém detalhes suficientes do desenho da coleta para sustentar afirmações de causa e efeito.