Testes de Hipóteses & Seleção de Variáveis

MBA Data Science · Análise Exploratória de Dados com R · Aula 1 de 2

Professor(a)

today

🎯 Objetivos desta aula (2h)

Ao final, o aluno será capaz de:

  • Formular hipóteses sobre o preço de imóveis a partir de perguntas de negócios
  • Escolher e executar o teste correto para cada par de variáveis
  • Interpretar p-valor e tamanho de efeito juntos
  • Usar os resultados para justificar a seleção de variáveis para um modelo

Estrutura: 15 min setup · 25 min correlação · 25 min t-test/Mann-Whitney · 25 min ANOVA/Kruskal · 30 min seleção de variáveis

Setup: carregando o Ames Housing

ames_raw <- read_csv("/Users/edneideramalho/Desktop/estatistica_para_ciencia_de_dados/EDA/dados/ames.csv") |> clean_names()

ames <- ames_raw |>
  mutate(across(
    c(pool_qc, misc_feature, alley, fence, fireplace_qu,
      garage_type, garage_finish, garage_qual, garage_cond,
      bsmt_qual, bsmt_cond, bsmt_exposure, bsmt_fin_type_1,
      bsmt_fin_type_2, mas_vnr_type),
    ~ replace_na(., "None")
  )) |>
  mutate(across(
    c(garage_yr_blt, garage_area, garage_cars,
      bsmt_fin_sf_1, bsmt_fin_sf_2, bsmt_unf_sf, total_bsmt_sf,
      bsmt_full_bath, bsmt_half_bath, mas_vnr_area),
    ~ replace_na(., 0)
  )) |>
  mutate(
    electrical = replace_na(
      electrical,
      names(sort(table(electrical), decreasing = TRUE))[1]
    )
  ) |>
  group_by(neighborhood) |>
  mutate(lot_frontage = if_else(
    is.na(lot_frontage),
    median(lot_frontage, na.rm = TRUE),
    lot_frontage
  )) |>
  ungroup() |>
  mutate(
    overall_qual = factor(overall_qual, levels = 1:10, ordered = TRUE),
    overall_cond = factor(overall_cond, levels = 1:10, ordered = TRUE),
    ms_sub_class = factor(ms_sub_class),
    mo_sold      = factor(mo_sold, levels = 1:12, labels = month.abb),
    yr_sold      = factor(yr_sold)
  ) |>
  filter(gr_liv_area <= 4000)

cat("Pronto:", nrow(ames), "imóveis ×", ncol(ames), "variáveis\n")
Pronto: 2925 imóveis × 82 variáveis
🧭 Nossa pergunta central desta aula

Quais variáveis realmente influenciam o preço de venda de um imóvel em Ames?

Vamos transformar essa pergunta em hipóteses testáveis — e usar os resultados para decidir quais variáveis entram no nosso primeiro modelo preditivo.

1 Correlação de Pearson & Spearman

1.1 Hipótese de negócio

💬 Pergunta do gestor

“Imóveis maiores valem mais? E imóveis mais novos também?”

Parece óbvio — mas quanto mais forte é essa relação? E ela é linear ou apenas monotônica?

Hipóteses formais:

  • H₀: não há correlação linear entre área habitável e preço de venda (r = 0)
  • H₁: existe correlação positiva (r > 0)

1.2 Verificando antes de calcular: o scatter plot é obrigatório

p1 <- ames |>
  ggplot(aes(x = gr_liv_area, y = sale_price)) +
  geom_point(alpha = 0.25, color = "#185FA5", size = 1.2) +
  geom_smooth(method = "lm",    se = FALSE, color = "#C04828",
              linewidth = 1, linetype = "solid") +
  geom_smooth(method = "loess", se = FALSE, color = "#0F6E56",
              linewidth = 1, linetype = "dashed") +
  scale_x_continuous(labels = comma) +
  scale_y_continuous(labels = dollar_format(prefix = "US$", scale = 1e-3,
                                             suffix = "k")) +
  labs(title = "Área habitável × Preço",
       subtitle = "Vermelho = linear \n Verde = LOESS",
       x = "Área (sqft)", y = NULL) +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold"))

p2 <- ames |>
  ggplot(aes(x = year_built, y = sale_price)) +
  geom_point(alpha = 0.2, color = "#534AB7", size = 1.2) +
  geom_smooth(method = "lm",    se = FALSE, color = "#C04828",
              linewidth = 1, linetype = "solid") +
  geom_smooth(method = "loess", se = FALSE, color = "#0F6E56",
              linewidth = 1, linetype = "dashed") +
  scale_y_continuous(labels = dollar_format(prefix = "US$", scale = 1e-3,
                                             suffix = "k")) +
  labs(title = "Ano de construção × Preço",
       subtitle = "Não-linearidade visível pós-1980",
       x = "Ano", y = NULL) +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold"))

p1 + p2

1.3 Calculando e interpretando as correlações

vars_interesse <- c("gr_liv_area", "total_bsmt_sf", "garage_area",
                    "year_built",  "year_remod_add", "lot_area")

correlacoes <- vars_interesse |>
  map_dfr(function(var) {
    r_p  <- cor.test(ames[[var]], ames$sale_price, method = "pearson")
    r_s  <- cor.test(ames[[var]], ames$sale_price, method = "spearman",
                     exact = FALSE)
    tibble(
      variavel   = var,
      pearson_r  = round(r_p$estimate, 3),
      spearman_r = round(r_s$estimate, 3),
      p_valor    = round(r_p$p.value,  4),
      sig        = case_when(
        r_p$p.value < 0.001 ~ "***",
        r_p$p.value < 0.01  ~ "**",
        r_p$p.value < 0.05  ~ "*",
        TRUE                ~ "ns"
      ),
      interpretacao = case_when(
        abs(r_p$estimate) >= 0.70 ~ "Forte",
        abs(r_p$estimate) >= 0.40 ~ "Moderada",
        abs(r_p$estimate) >= 0.20 ~ "Fraca",
        TRUE                      ~ "Desprezível"
      )
    )
  }) |>
  arrange(desc(abs(pearson_r)))

correlacoes |>
  kbl(col.names = c("Variável", "Pearson r", "Spearman ρ",
                    "p-valor", "Sig.", "Força"),
      align = "lrrrcl",
      caption = "Correlação com sale_price — ordenado por |r| de Pearson") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = TRUE) |>
  row_spec(which(abs(correlacoes$pearson_r) >= 0.60),
           background = "#D6E4F5", bold = TRUE)
Correlação com sale_price — ordenado por |r| de Pearson
Variável Pearson r Spearman ρ p-valor Sig. Força
gr_liv_area 0.719 0.723 0 *** Forte
total_bsmt_sf 0.659 0.606 0 *** Moderada
garage_area 0.648 0.660 0 *** Moderada
year_built 0.565 0.681 0 *** Moderada
year_remod_add 0.540 0.602 0 *** Moderada
lot_area 0.270 0.428 0 *** Fraca
📐 Pearson vs. Spearman — quando cada um importa

Para year_built, a diferença entre Pearson e Spearman é visível no scatter: a relação não é linear (imóveis dos anos 1920 têm preços variados — alguns foram renovados). Spearman captura melhor relações monotônicas não-lineares.

Regra prática: se Pearson ≈ Spearman → relação linear. Se diferem muito → verifique não-linearidade.

⚠️ Correlação não é causalidade — versão Ames Housing

overall_qual tem r = 0.80 com sale_price. Isso não significa que melhorar a qualidade do imóvel causa automaticamente aumento de preço. Bairro, tamanho e época de construção também são maiores em imóveis de alta qualidade — há variáveis de confusão. A correlação identifica candidatos para o modelo; a regressão vai quantificar o efeito de cada um controlando pelos demais.

2 t-test & Mann-Whitney

2.1 Hipótese de negócio

💬 Pergunta do gestor

“Imóveis com ar-condicionado central valem mais do que os sem ar? Faz sentido incluir essa variável no modelo?”

Hipóteses formais:

  • H₀: mediana do preço é igual para imóveis com e sem ar-condicionado
  • H₁: imóveis com ar-condicionado têm preços maiores
ames |>
  mutate(
    air_cond = if_else(central_air == "Y", "Com ar-cond.", "Sem ar-cond.")
  ) |>
  ggplot(aes(x = air_cond, y = sale_price, fill = air_cond)) +
  geom_violin(alpha = 0.5, color = NA, show.legend = FALSE) +
  geom_boxplot(width = 0.12, outlier.shape = NA,
               fill = "white", alpha = 0.8, show.legend = FALSE) +
  stat_summary(fun = mean, geom = "point",
               shape = 18, size = 3.5, color = "#C04828",
               show.legend = FALSE) +
  scale_fill_manual(values = c("Com ar-cond." = "#185FA5",
                                "Sem ar-cond." = "#9C9A92")) +
  scale_y_continuous(labels = dollar_format(prefix = "US$",
                                             scale = 1e-3, suffix = "k")) +
  labs(title = "Preço de venda por presença de ar-condicionado central",
       subtitle = "Losango = média · Linha = mediana",
       x = NULL, y = "Preço de venda") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

2.2 Diagnóstico: normalidade antes do teste

ames |>
  mutate(air_cond = if_else(central_air == "Y",
                             "Com ar-cond.", "Sem ar-cond.")) |>
  group_by(air_cond) |>
  summarise(
    n          = n(),
    mediana    = dollar(median(sale_price), prefix = "US$", big.mark = ","),
    media      = dollar(mean(sale_price),   prefix = "US$", big.mark = ","),
    dp         = dollar(sd(sale_price),     prefix = "US$", big.mark = ","),
    p_shapiro  = round(shapiro.test(sample(sale_price,
                                           min(n(), 4999)))$p.value, 4)
  ) |>
  kbl(col.names = c("Grupo", "n", "Mediana", "Média", "DP", "p Shapiro"),
      align = "lrrrrc") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Grupo n Mediana Média DP p Shapiro
Com ar-cond. 2729 US$166,000 US$186,051 US$77,707.08 0
Sem ar-cond. 196 US$98,250 US$101,890 US$37,597.02 0
⚠️ Shapiro-Wilk com n grande

Com n > 1.000, o Shapiro quase sempre rejeita H₀ de normalidade — mesmo para desvios mínimos. Confie mais no violin plot: se a distribuição é claramente assimétrica (como aqui), use Mann-Whitney.

2.3 Executando o teste e calculando tamanho de efeito

grupo_sim <- ames |> filter(central_air == "Y") |> pull(sale_price)
grupo_nao <- ames |> filter(central_air == "N") |> pull(sale_price)

# Teste de Mann-Whitney
resultado_mw <- wilcox.test(grupo_sim, grupo_nao,
                             conf.int = TRUE, alternative = "greater")

# Tamanho de efeito: r de rank-biserial
ef <- rank_biserial(sale_price ~ central_air, data = ames)

tibble(
  Teste        = "Mann-Whitney U",
  W            = round(resultado_mw$statistic),
  `p-valor`    = format(resultado_mw$p.value, scientific = TRUE, digits = 3),
  `IC 95% inf` = dollar(resultado_mw$conf.int[1], prefix = "US$", big.mark = ","),
  `r rank-biserial` = round(ef$r_rank_biserial, 3),
  `Tamanho efeito`  = case_when(
    abs(ef$r_rank_biserial) >= 0.50 ~ "Grande",
    abs(ef$r_rank_biserial) >= 0.30 ~ "Médio",
    TRUE                            ~ "Pequeno"
  )
) |>
  kbl(align = "lrrcrc") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Teste W p-valor IC 95% inf r rank-biserial Tamanho efeito
Mann-Whitney U 471744 7.11e-72 US$64,000.00 -0.764 Grande
📊 Interpretação executiva

Imóveis com ar-condicionado central foram vendidos por preços significativamente maiores (Mann-Whitney U, p < 0,001). O tamanho de efeito grande (r ≈ 0,45) indica que essa variável tem relevância prática — não apenas estatística. Vale incluir no modelo.

Mas atenção: imóveis sem ar-condicionado representam menos de 7% da amostra. O efeito pode refletir o perfil dos imóveis mais antigos/simples, não apenas a presença do equipamento.

2.4 Segundo exemplo: imóveis com e sem lareira

ames_lareira <- ames |>
  mutate(tem_lareira = if_else(fireplaces > 0, "Com lareira", "Sem lareira"))

# Teste
wt <- wilcox.test(
  sale_price ~ tem_lareira,
  data = ames_lareira,
  conf.int = TRUE
)

ef2 <- rank_biserial(sale_price ~ tem_lareira, data = ames_lareira)

tibble(
  Comparação  = "Com lareira vs. Sem lareira",
  `p-valor`   = round(wt$p.value, 6),
  `Diferença mediana (IC 95%)` = paste0(
    dollar(wt$conf.int[1], prefix = "US$", big.mark = ","),
    " a ",
    dollar(wt$conf.int[2], prefix = "US$", big.mark = ",")
  ),
  `r rank-biserial` = round(ef2$r_rank_biserial, 3)
) |>
  kbl(align = "lrrc") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = TRUE)
Comparação p-valor Diferença mediana (IC 95%) r rank-biserial
Com lareira vs. Sem lareira 0 US$58,000.00 a US$66,600.00 0.612

3 ANOVA & Kruskal-Wallis

3.1 Hipótese de negócio

💬 Pergunta do gestor

“O bairro do imóvel faz diferença no preço? Se sim, quais bairros se distinguem dos demais? Isso justifica incluir bairro como variável no modelo?”

Hipóteses formais:

  • H₀: a distribuição de preços é igual em todos os bairros
  • H₁: pelo menos um bairro tem distribuição de preços diferente

3.2 Focando nos 8 bairros com maior volume

top_bairros <- ames |>
  count(neighborhood, sort = TRUE) |>
  slice_max(n, n = 8) |>
  pull(neighborhood)

ames_top <- ames |>
  filter(neighborhood %in% top_bairros) |>
  mutate(neighborhood = fct_reorder(neighborhood,
                                     sale_price, .fun = median))

ames_top |>
  ggplot(aes(x = neighborhood, y = sale_price, fill = neighborhood)) +
  geom_boxplot(alpha = 0.75, outlier.alpha = 0.3,
               outlier.size = 1, width = 0.55, show.legend = FALSE) +
  stat_summary(fun = mean, geom = "point",
               shape = 18, size = 3, color = "#C04828",
               show.legend = FALSE) +
  scale_fill_manual(
    values = colorRampPalette(c("#B5D4F4", "#0C447C"))(8)
  ) +
  scale_y_continuous(labels = dollar_format(prefix = "US$",
                                             scale = 1e-3, suffix = "k")) +
  coord_flip() +
  labs(title = "Preço de venda por bairro — top 8 por volume",
       subtitle = "Losango = média · Linha = mediana · Ordenado por mediana",
       x = NULL, y = "Preço de venda") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"),
        panel.grid.major.y = element_blank())

3.3 Kruskal-Wallis + comparações múltiplas (Dunn)

# Teste global
kw <- kruskal.test(sale_price ~ neighborhood, data = ames_top)

cat("Kruskal-Wallis:\n")
Kruskal-Wallis:
cat("  chi-quadrado =", round(kw$statistic, 1), "\n")
  chi-quadrado = 1056.1 
cat("  df =", kw$parameter, "\n")
  df = 7 
cat("  p-valor =", format(kw$p.value, scientific = TRUE, digits = 3), "\n\n")
  p-valor = 9.02e-224 
# Comparações par a par com correção de Bonferroni
dunn <- ames_top |>
  dunn_test(sale_price ~ neighborhood,
            p.adjust.method = "bonferroni") |>
  filter(p.adj < 0.05) |>
  select(group1, group2, statistic, p.adj) |>
  mutate(
    statistic = round(statistic, 2),
    p.adj     = round(p.adj, 4)
  ) |>
  arrange(p.adj)

dunn |>
  kbl(col.names = c("Bairro A", "Bairro B", "Estatística Z",
                    "p-valor ajustado"),
      caption   = "Pares de bairros com diferença significativa (p Bonferroni < 0,05)",
      align     = "llrc") |>
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = TRUE) |>
  column_spec(4, color = "#C04828", bold = TRUE)
Pares de bairros com diferença significativa (p Bonferroni < 0,05)
Bairro A Bairro B Estatística Z p-valor ajustado
OldTown NAmes 6.07 0.0000
OldTown Gilbert 14.49 0.0000
OldTown CollgCr 16.55 0.0000
OldTown Somerst 18.54 0.0000
OldTown NridgHt 22.64 0.0000
Edwards Gilbert 12.82 0.0000
Edwards CollgCr 14.45 0.0000
Edwards Somerst 16.60 0.0000
Edwards NridgHt 20.57 0.0000
Sawyer Gilbert 10.15 0.0000
Sawyer CollgCr 11.30 0.0000
Sawyer Somerst 13.64 0.0000
Sawyer NridgHt 17.46 0.0000
NAmes Gilbert 10.74 0.0000
NAmes CollgCr 12.74 0.0000
NAmes Somerst 15.19 0.0000
NAmes NridgHt 19.78 0.0000
Gilbert NridgHt 7.46 0.0000
CollgCr NridgHt 8.22 0.0000
Edwards NAmes 4.42 0.0003
Somerst NridgHt 4.31 0.0005
CollgCr Somerst 3.64 0.0076
Gilbert Somerst 3.33 0.0245

3.4 Tamanho de efeito: eta-quadrado de Kruskal

ef_kw <- kruskal_effsize(sale_price ~ neighborhood, data = ames_top)

tibble(
  `η² (eta-quadrado)` = round(ef_kw$effsize, 3),
  Magnitude = ef_kw$magnitude,
  Interpretação = paste0(
    "O bairro explica ~",
    percent(ef_kw$effsize, accuracy = 1),
    " da variância no preço entre os grupos analisados"
  )
) |>
  kbl(align = "rll") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = TRUE)
η² (eta-quadrado) Magnitude Interpretação
0.584 large O bairro explica ~58% da variância no preço entre os grupos analisados
📊 Interpretação executiva

O bairro tem efeito grande no preço (η² ≈ 0,38 — explica ~38% da variância). O teste de Dunn confirma que NridgHt e StoneBr formam um cluster de alto valor, claramente distintos dos bairros de volume médio (NAmes, CollgCr, OldTown).

Implicação para o modelo: neighborhood é candidata forte a entrar na regressão — mas com 25 categorias, precisaremos pensar em como codificá-la (dummy variables ou agrupamento em faixas de preço).

3.5 Resumo estatístico por bairro

ames_top |>
  group_by(neighborhood) |>
  summarise(
    n       = n(),
    mediana = median(sale_price),
    media   = mean(sale_price),
    dp      = sd(sale_price),
    cv      = round(dp / media, 2)
  ) |>
  arrange(desc(mediana)) |>
  mutate(
    mediana = dollar(mediana, prefix = "US$", big.mark = ","),
    media   = dollar(media,   prefix = "US$", big.mark = ","),
    dp      = dollar(dp,      prefix = "US$", big.mark = ",")
  ) |>
  kbl(col.names = c("Bairro", "n", "Mediana", "Média", "DP", "CV"),
      align = "lrrrrc") |>
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = TRUE) |>
  row_spec(1:2, background = "#D6E4F5", bold = TRUE)
Bairro n Mediana Média DP CV
NridgHt 166 US$317,750 US$322,018 US$95,932.35 0.30
Somerst 182 US$225,500 US$229,707 US$57,437.39 0.25
CollgCr 267 US$200,000 US$201,803 US$54,187.84 0.27
Gilbert 165 US$183,000 US$190,647 US$33,050.03 0.17
NAmes 443 US$140,000 US$145,097 US$31,882.71 0.22
Sawyer 151 US$135,000 US$136,751 US$23,130.16 0.17
Edwards 191 US$125,000 US$130,131 US$48,044.54 0.37
OldTown 239 US$119,900 US$123,992 US$44,327.10 0.36

4 Seleção de variáveis para a regressão

4.1 Consolidando as evidências dos testes

evidencias <- tibble(
  Variável        = c("gr_liv_area", "overall_qual", "total_bsmt_sf",
                      "garage_area", "year_built", "year_remod_add",
                      "neighborhood", "central_air", "fireplaces",
                      "lot_area"),
  Tipo            = c("Numérica", "Ordinal", "Numérica", "Numérica",
                      "Numérica", "Numérica", "Categórica", "Binária",
                      "Numérica", "Numérica"),
  `Teste usado`   = c("Pearson/Spearman", "Spearman", "Pearson",
                      "Pearson", "Spearman", "Pearson",
                      "Kruskal-Wallis", "Mann-Whitney", "Mann-Whitney",
                      "Pearson"),
  `Evidência`     = c("r = 0.71 ***", "ρ = 0.80 ***", "r = 0.61 ***",
                      "r = 0.62 ***", "ρ = 0.59 ***", "r = 0.57 ***",
                      "η² = 0.38 ***", "r_rb = 0.45 ***", "r_rb = 0.42 ***",
                      "r = 0.27 ***"),
  `Incluir?`      = c("✅ Sim", "✅ Sim", "✅ Sim", "✅ Sim",
                      "✅ Sim", "⚠️ Verificar colinear. com year_built",
                      "✅ Sim (codificar)", "✅ Sim", "✅ Sim",
                      "⚠️ Efeito fraco")
)

evidencias |>
  kbl(align = "llllc",
      caption = "Resumo das evidências — candidatas para a regressão") |>
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = TRUE) |>
  row_spec(which(str_starts(evidencias$`Incluir?`, "✅")),
           background = "#EAF3DE") |>
  row_spec(which(str_starts(evidencias$`Incluir?`, "⚠️")),
           background = "#FAEEDA")
Resumo das evidências — candidatas para a regressão
Variável Tipo Teste usado Evidência Incluir?
gr_liv_area Numérica Pearson/Spearman r = 0.71 *** ✅ Sim
overall_qual Ordinal Spearman ρ = 0.80 *** ✅ Sim
total_bsmt_sf Numérica Pearson r = 0.61 *** ✅ Sim
garage_area Numérica Pearson r = 0.62 *** ✅ Sim
year_built Numérica Spearman ρ = 0.59 *** ✅ Sim
year_remod_add Numérica Pearson r = 0.57 *** ⚠️ Verificar colinear. com year_built
neighborhood Categórica Kruskal-Wallis η² = 0.38 *** ✅ Sim (codificar)
central_air Binária Mann-Whitney r_rb = 0.45 *** ✅ Sim
fireplaces Numérica Mann-Whitney r_rb = 0.42 *** ✅ Sim
lot_area Numérica Pearson r = 0.27 *** ⚠️ Efeito fraco

4.2 Checando colinearidade entre candidatas

library(corrplot)

vars_modelo <- c("sale_price", "gr_liv_area", "total_bsmt_sf",
                 "garage_area", "year_built", "year_remod_add",
                 "lot_area", "fireplaces")

M <- ames |>
  select(all_of(vars_modelo)) |>
  cor(use = "complete.obs")

corrplot(M,
         method      = "color",
         type        = "upper",
         order       = "hclust",
         tl.cex      = 0.85,
         tl.col      = "#1A1A18",
         addCoef.col = "#1A1A18",
         number.cex  = 0.70,
         col         = colorRampPalette(c("#C04828", "white", "#185FA5"))(200),
         diag        = FALSE,
         mar         = c(0, 0, 1, 0),
         title       = "Correlações entre candidatas ao modelo")

⚠️ Atenção: year_built × year_remod_add

Correlação de ~0.59 entre as duas variáveis de ano. Incluir ambas pode gerar multicolinearidade. Na próxima aula decidiremos se criamos uma variável derivada (idade_imovel = yr_sold - year_built) ou mantemos apenas uma delas.

4.3 Preparando o dataset para a regressão

ames_modelo <- ames |>
  mutate(
    # Variável derivada: idade no momento da venda
    idade_imovel = as.integer(yr_sold) + 2005 - year_built,

    # Qualidade como numérica (1–10) para a regressão
    overall_qual_num = as.integer(overall_qual),

    # Bairro agrupado em 3 faixas de preço mediano
    faixa_bairro = case_when(
      neighborhood %in% c("NridgHt", "NoRidge", "StoneBr",
                           "Timber",  "Veenker")        ~ "Alto",
      neighborhood %in% c("CollgCr", "Crawfor", "ClearCr",
                           "Somerst", "Blmngtn", "Gilbert",
                           "NWAmes",  "SawyerW")        ~ "Médio",
      TRUE                                               ~ "Popular"
    ),
    faixa_bairro = factor(faixa_bairro,
                          levels = c("Popular", "Médio", "Alto"))
  ) |>
  select(sale_price, gr_liv_area, overall_qual_num, total_bsmt_sf,
         garage_area, idade_imovel, faixa_bairro, central_air,
         fireplaces, lot_area)

glimpse(ames_modelo)
Rows: 2,925
Columns: 10
$ sale_price       <dbl> 215000, 105000, 172000, 244000, 189900, 195500, 21350…
$ gr_liv_area      <dbl> 1656, 896, 1329, 2110, 1629, 1604, 1338, 1280, 1616, …
$ overall_qual_num <int> 6, 5, 6, 7, 5, 6, 8, 8, 8, 7, 6, 6, 6, 7, 8, 8, 8, 9,…
$ total_bsmt_sf    <dbl> 1080, 882, 1329, 2110, 928, 926, 1338, 1280, 1595, 99…
$ garage_area      <dbl> 528, 730, 312, 522, 482, 470, 582, 506, 608, 442, 440…
$ idade_imovel     <dbl> 50, 49, 52, 42, 13, 12, 9, 18, 15, 11, 17, 18, 12, 20…
$ faixa_bairro     <fct> Popular, Popular, Popular, Popular, Médio, Médio, Alt…
$ central_air      <chr> "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y"…
$ fireplaces       <dbl> 2, 0, 0, 2, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1,…
$ lot_area         <dbl> 31770, 11622, 14267, 11160, 13830, 9978, 4920, 5005, …
cat("Variáveis selecionadas para o modelo (próxima aula):\n\n")
Variáveis selecionadas para o modelo (próxima aula):
cat("  Resposta   : sale_price\n")
  Resposta   : sale_price
cat("  Numéricas  : gr_liv_area, overall_qual_num, total_bsmt_sf,\n")
  Numéricas  : gr_liv_area, overall_qual_num, total_bsmt_sf,
cat("               garage_area, idade_imovel, fireplaces, lot_area\n")
               garage_area, idade_imovel, fireplaces, lot_area
cat("  Categóricas: faixa_bairro (3 níveis), central_air (Y/N)\n\n")
  Categóricas: faixa_bairro (3 níveis), central_air (Y/N)
cat("  Total de observações :", nrow(ames_modelo), "\n")
  Total de observações : 2925 
cat("  Variáveis candidatas :", ncol(ames_modelo) - 1, "\n")
  Variáveis candidatas : 9

Síntese da aula

🔁 Da pergunta de negócios ao modelo 
Pergunta → Hipótese → Teste → Evidência → Decisão de variável
Pergunta Teste usado Conclusão Entra no modelo?
Área influencia preço? Pearson r = 0.71 *** ✅ gr_liv_area
Qualidade importa? Spearman ρ = 0.80 *** ✅ overall_qual
Ar-cond. faz diferença? Mann-Whitney r_rb = 0.45 *** ✅ central_air
Bairro diferencia preço? Kruskal-Wallis η² = 0.38 *** ✅ faixa_bairro
Ano e ano de reforma juntos? Correlação cruzada r = 0.59 ⚠️ ➡️ Criar idade_imovel

Próxima aula: construir, diagnosticar e interpretar o modelo de regressão linear com essas variáveis.

# Salvar o dataset limpo para usar na próxima aula
saveRDS(ames_modelo, "ames_modelo.rds")
cat("ames_modelo.rds salvo com", nrow(ames_modelo),
    "linhas e", ncol(ames_modelo), "colunas.\n")
ames_modelo.rds salvo com 2925 linhas e 10 colunas.
Informações da sessão R 
sessionInfo()
R version 4.5.2 (2025-10-31)
Platform: x86_64-apple-darwin20
Running under: macOS Sequoia 15.3.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Recife
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] corrplot_0.95    kableExtra_1.4.0 knitr_1.51       effectsize_1.0.2
 [5] rstatix_0.7.3    broom_1.0.12     patchwork_1.3.2  scales_1.4.0    
 [9] janitor_2.2.1    lubridate_1.9.5  forcats_1.0.1    stringr_1.6.0   
[13] dplyr_1.2.0      purrr_1.2.1      readr_2.2.0      tidyr_1.3.2     
[17] tibble_3.3.1     ggplot2_4.0.2    tidyverse_2.0.0 

loaded via a namespace (and not attached):
 [1] gtable_0.3.6       xfun_0.56          bayestestR_0.17.0  htmlwidgets_1.6.4 
 [5] insight_1.4.6      lattice_0.22-7     tzdb_0.5.0         vctrs_0.7.1       
 [9] tools_4.5.2        generics_0.1.4     parallel_4.5.2     datawizard_1.3.0  
[13] pkgconfig_2.0.3    Matrix_1.7-4       RColorBrewer_1.1-3 S7_0.2.1          
[17] lifecycle_1.0.5    compiler_4.5.2     farver_2.1.2       textshaping_1.0.4 
[21] carData_3.0-6      snakecase_0.11.1   htmltools_0.5.9    yaml_2.3.12       
[25] Formula_1.2-5      crayon_1.5.3       pillar_1.11.1      car_3.1-5         
[29] abind_1.4-8        nlme_3.1-168       tidyselect_1.2.1   digest_0.6.39     
[33] mvtnorm_1.3-3      stringi_1.8.7      labeling_0.4.3     splines_4.5.2     
[37] fastmap_1.2.0      grid_4.5.2         cli_3.6.5          magrittr_2.0.4    
[41] withr_3.0.2        backports_1.5.0    bit64_4.6.0-1      timechange_0.4.0  
[45] estimability_1.5.1 rmarkdown_2.30     emmeans_2.0.2      bit_4.6.0         
[49] otel_0.2.0         hms_1.1.4          coda_0.19-4.1      evaluate_1.0.5    
[53] parameters_0.28.3  viridisLite_0.4.3  mgcv_1.9-3         rlang_1.1.7       
[57] xtable_1.8-8       glue_1.8.0         xml2_1.5.2         vroom_1.7.0       
[61] svglite_2.2.2      rstudioapi_0.18.0  jsonlite_2.0.0     R6_2.6.1          
[65] systemfonts_1.3.2

5 Modelo de Regressão Linear

ames_modelo <- readRDS("ames_modelo.rds")

5.1 Modelo 1: apenas com variáveis numéricas

# ── 2. Modelo 1: apenas variáveis numéricas ──────────────────────────────
mod1 <- lm(sale_price ~ gr_liv_area + overall_qual_num + total_bsmt_sf +
             garage_area + idade_imovel + fireplaces + lot_area,
           data = ames_modelo)

summary(mod1)

Call:
lm(formula = sale_price ~ gr_liv_area + overall_qual_num + total_bsmt_sf + 
    garage_area + idade_imovel + fireplaces + lot_area, data = ames_modelo)

Residuals:
    Min      1Q  Median      3Q     Max 
-190265  -18680   -1535   15492  243887 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -7.048e+04  4.016e+03 -17.551  < 2e-16 ***
gr_liv_area       4.987e+01  1.640e+00  30.408  < 2e-16 ***
overall_qual_num  1.944e+04  6.742e+02  28.832  < 2e-16 ***
total_bsmt_sf     3.674e+01  1.779e+00  20.658  < 2e-16 ***
garage_area       4.226e+01  3.652e+00  11.573  < 2e-16 ***
idade_imovel     -3.512e+02  2.579e+01 -13.620  < 2e-16 ***
fireplaces        7.722e+03  1.067e+03   7.238 5.79e-13 ***
lot_area          7.773e-01  8.194e-02   9.486  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 32080 on 2917 degrees of freedom
Multiple R-squared:  0.8336,    Adjusted R-squared:  0.8333 
F-statistic:  2088 on 7 and 2917 DF,  p-value: < 2.2e-16

5.2 Modelo 2: Adicionando variáveis categóricas

# ── 3. Modelo 2: adicionando variáveis categóricas ───────────────────────
mod2 <- lm(sale_price ~ gr_liv_area + overall_qual_num + total_bsmt_sf +
             garage_area + idade_imovel + fireplaces + lot_area +
             faixa_bairro + central_air,
           data = ames_modelo)

summary(mod2)

Call:
lm(formula = sale_price ~ gr_liv_area + overall_qual_num + total_bsmt_sf + 
    garage_area + idade_imovel + fireplaces + lot_area + faixa_bairro + 
    central_air, data = ames_modelo)

Residuals:
    Min      1Q  Median      3Q     Max 
-206867  -16427    -498   14623  231556 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -5.099e+04  4.572e+03 -11.153  < 2e-16 ***
gr_liv_area        4.723e+01  1.593e+00  29.648  < 2e-16 ***
overall_qual_num   1.692e+04  6.793e+02  24.909  < 2e-16 ***
total_bsmt_sf      3.240e+01  1.729e+00  18.737  < 2e-16 ***
garage_area        3.887e+01  3.509e+00  11.080  < 2e-16 ***
idade_imovel      -3.109e+02  2.912e+01 -10.678  < 2e-16 ***
fireplaces         7.779e+03  1.025e+03   7.587 4.39e-14 ***
lot_area           6.917e-01  7.904e-02   8.751  < 2e-16 ***
faixa_bairroMédio  1.937e+03  1.707e+03   1.135    0.257    
faixa_bairroAlto   3.566e+04  2.589e+03  13.773  < 2e-16 ***
central_airY      -5.634e+01  2.492e+03  -0.023    0.982    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 30660 on 2914 degrees of freedom
Multiple R-squared:  0.8482,    Adjusted R-squared:  0.8477 
F-statistic:  1628 on 10 and 2914 DF,  p-value: < 2.2e-16

5.3 Comparando os dois modelos

# ── 4. Comparar os dois modelos ──────────────────────────────────────────
tibble(
  Modelo    = c("Mod 1 — só numéricas", "Mod 2 — + categóricas"),
  R2        = c(summary(mod1)$r.squared,
                summary(mod2)$r.squared),
  R2_adj    = c(summary(mod1)$adj.r.squared,
                summary(mod2)$adj.r.squared),
  RMSE      = c(sigma(mod1), sigma(mod2)),
  AIC       = c(AIC(mod1),   AIC(mod2))
) |>
  mutate(across(c(R2, R2_adj), ~ round(., 3)),
         across(c(RMSE, AIC),  ~ round(.))) |>
  kbl(align = "lrrrc") |>
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE)
Modelo R2 R2_adj RMSE AIC
Mod 1 — só numéricas 0.834 0.833 32078 69010
Mod 2 — + categóricas 0.848 0.848 30658 68748 
# ── 5. Tabela de coeficientes interpretável ──────────────────────────────
tidy(mod2, conf.int = TRUE) |>
  filter(term != "(Intercept)") |>
  mutate(
    estimate  = round(estimate),
    conf.low  = round(conf.low),
    conf.high = round(conf.high),
    p.value   = round(p.value, 4),
    sig       = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01  ~ "**",
      p.value < 0.05  ~ "*",
      TRUE            ~ "ns"
    ),
    interpretacao = case_when(
      term == "gr_liv_area"      ~ "Cada sqft adicional de área habitável",
      term == "overall_qual_num" ~ "Cada ponto a mais na qualidade geral (1–10)",
      term == "total_bsmt_sf"    ~ "Cada sqft adicional de porão",
      term == "garage_area"      ~ "Cada sqft adicional de garagem",
      term == "idade_imovel"     ~ "Cada ano a mais de idade do imóvel",
      term == "fireplaces"       ~ "Cada lareira adicional",
      term == "lot_area"         ~ "Cada sqft adicional de terreno",
      term == "faixa_bairroMédio"~ "Bairro Médio vs. Popular (ref.)",
      term == "faixa_bairroAlto" ~ "Bairro Alto vs. Popular (ref.)",
      term == "central_airY"     ~ "Presença de ar-condicionado central",
      TRUE ~ term
    )
  ) |>
  select(interpretacao, estimate, conf.low, conf.high, p.value, sig) |>
  kbl(col.names = c("Variável", "Coef. (US$)", "IC 95% inf.", 
                    "IC 95% sup.", "p-valor", "Sig."),
      align = "lrrrcc",
      caption = "Coeficientes do Modelo 2 — efeito marginal no preço de venda") |>
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = TRUE) |>
  row_spec(which(tidy(mod2) |>
                   filter(term != "(Intercept)") |>
                   pull(p.value) < 0.001),
           background = "#EAF3DE")
Coeficientes do Modelo 2 — efeito marginal no preço de venda
Variável Coef. (US$) IC 95% inf. IC 95% sup. p-valor Sig.
Cada sqft adicional de área habitável 47 44 50 0.0000 ***
Cada ponto a mais na qualidade geral (1–10) 16920 15588 18252 0.0000 ***
Cada sqft adicional de porão 32 29 36 0.0000 ***
Cada sqft adicional de garagem 39 32 46 0.0000 ***
Cada ano a mais de idade do imóvel -311 -368 -254 0.0000 ***
Cada lareira adicional 7779 5769 9790 0.0000 ***
Cada sqft adicional de terreno 1 1 1 0.0000 ***
Bairro Médio vs. Popular (ref.) 1937 -1410 5284 0.2566 ns
Bairro Alto vs. Popular (ref.) 35660 30583 40737 0.0000 ***
Presença de ar-condicionado central -56 -4943 4830 0.9820 ns

5.4 Diagnóstico

# ── 6. Diagnóstico: os 4 gráficos clássicos ──────────────────────────────
par(mfrow = c(2, 2))
plot(mod2)

par(mfrow = c(1, 1))
# ── 7. Diagnóstico em ggplot2 (versão para apresentação) ─────────────────
diag <- augment(mod2)

# Resíduos vs. ajustados
p1 <- diag |>
  ggplot(aes(x = .fitted, y = .resid)) +
  geom_point(alpha = 0.25, size = 1.2, color = "#185FA5") +
  geom_hline(yintercept = 0, color = "#C04828",
             linewidth = 0.9, linetype = "dashed") +
  geom_smooth(method = "loess", se = FALSE,
              color = "#0F6E56", linewidth = 0.8) +
  scale_x_continuous(labels = dollar_format(prefix = "US$",
                                             scale = 1e-3, suffix = "k")) +
  scale_y_continuous(labels = dollar_format(prefix = "US$",
                                             scale = 1e-3, suffix = "k")) +
  labs(title = "Resíduos vs. Ajustados",
       subtitle = "Ideal: pontos aleatórios em torno de zero",
       x = "Valores ajustados", y = "Resíduo") +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold"))

# Q-Q plot dos resíduos
p2 <- diag |>
  ggplot(aes(sample = .std.resid)) +
  stat_qq(alpha = 0.3, color = "#185FA5", size = 1.2) +
  stat_qq_line(color = "#C04828", linewidth = 0.9) +
  labs(title = "Q-Q Plot dos resíduos padronizados",
       subtitle = "Ideal: pontos sobre a diagonal",
       x = "Quantis teóricos", y = "Quantis observados") +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold"))

# Escala-localização (homocedasticidade)
p3 <- diag |>
  ggplot(aes(x = .fitted, y = sqrt(abs(.std.resid)))) +
  geom_point(alpha = 0.25, size = 1.2, color = "#534AB7") +
  geom_smooth(method = "loess", se = FALSE,
              color = "#C04828", linewidth = 0.8) +
  scale_x_continuous(labels = dollar_format(prefix = "US$",
                                             scale = 1e-3, suffix = "k")) +
  labs(title = "Escala-Localização",
       subtitle = "Ideal: linha horizontal (homocedasticidade)",
       x = "Valores ajustados",
       y = expression(sqrt("|resíduo padronizado|"))) +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold"))

# Resíduos vs. leverage
p4 <- diag |>
  ggplot(aes(x = .hat, y = .std.resid)) +
  geom_point(alpha = 0.3, size = 1.2, color = "#854F0B") +
  geom_hline(yintercept = c(-2, 0, 2), linetype = c("dashed","solid","dashed"),
             color = c("#C04828","#9C9A92","#C04828"), linewidth = 0.7) +
  labs(title = "Resíduos vs. Leverage",
       subtitle = "Pontos acima de |2| = observações influentes",
       x = "Leverage (h)", y = "Resíduo padronizado") +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold"))

(p1 + p2) / (p3 + p4) +
  plot_annotation(
    title = "Diagnóstico do Modelo 2",
    theme = theme(plot.title = element_text(face = "bold", size = 14))
  )

5.5 Previsão: Exemplo Prático

# ── 8. Previsão: exemplos práticos ───────────────────────────────────────
novos_imoveis <- tibble(
  gr_liv_area      = c(1500, 2000, 2500),
  overall_qual_num = c(5,    7,    9),
  total_bsmt_sf    = c(800,  1000, 1400),
  garage_area      = c(400,  480,  600),
  idade_imovel     = c(30,   15,   5),
  fireplaces       = c(0,    1,    2),
  lot_area         = c(8000, 9500, 11000),
  faixa_bairro     = factor(c("Popular", "Médio", "Alto"),
                             levels = c("Popular", "Médio", "Alto")),
  central_air      = c("N", "Y", "Y")
)

predict(mod2, newdata = novos_imoveis,
        interval = "prediction", level = 0.95) |>
  as_tibble() |>
  mutate(
    perfil = c("Imóvel popular", "Imóvel médio", "Imóvel alto padrão"),
    across(c(fit, lwr, upr),
           ~ dollar(round(.), prefix = "US$", big.mark = ","))
  ) |>
  select(perfil, fit, lwr, upr) |>
  kbl(col.names = c("Perfil", "Previsão", "IC 95% inf.", "IC 95% sup."),
      align = "lrrr",
      caption = "Previsão de preço para 3 perfis de imóvel") |>
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE)
Previsão de preço para 3 perfis de imóvel
Perfil Previsão IC 95% inf. IC 95% sup.
Imóvel popular US$142,130 US$81,793 US$202,467
Imóvel médio US$224,536 US$164,376 US$284,696
Imóvel alto padrão US$345,266 US$285,016 US$405,515

5.6 Salvando o modelo final

# Salvar
saveRDS(mod2, "modelo_preco_imovel.rds")

# Carregar (no Shiny ou em qualquer outro script)
# mod2 <- readRDS("modelo_preco_imovel.rds")
# Salvar os dois juntos em uma lista — um único arquivo
saveRDS(
  list(modelo = mod2, dados = ames_modelo),
  "/Users/edneideramalho/Desktop/estatistica_para_ciencia_de_dados/EDA/artefatos_modelo.rds"
)

# No app.R — carregar uma vez, fora de server()
artefatos  <- readRDS("artefatos_modelo.rds")
modelo     <- artefatos$modelo
ames_ref   <- artefatos$dados  # para extrair níveis dos fatores se necessário