Análise de Dados e Regressão Linear — Qualidade do Ar

Autor

Lucas Levy Peruchi

Data de Publicação

17 de junho de 2026

Código publicado no RPubs: https://rpubs.com/llevyperuchi/ozonio

Código

library(tidyverse)   # dplyr + ggplot2
library(gtsummary)   # tbl_summary

1. A base de dados escolhida

Escolhi a base airquality, que já vem no R. Ela tem medições diárias de qualidade do ar em Nova York entre maio e setembro de 1973. Escolhi ela por dois motivos práticos: já tem dados faltantes de verdade (nas colunas Ozone e Solar.R), que é uma exigência do trabalho, e tem variáveis numéricas que fazem sentido relacionar (clima x poluição).

A base original tem só Month e Day como possíveis categóricas, e elas são fracas (números). Pra cumprir a exigência de pelo menos 2 variáveis categóricas, eu derivei duas:

Mes: o mês como fator nomeado (Maio…Setembro);
Faixa_Temp: faixa de temperatura (Fria / Amena / Quente), cortando a coluna Temp.

Minha variável de interesse é o Ozone (nível de ozônio), que é numérica.

O que eu espero encontrar: que o ozônio sobe quando a temperatura sobe e cai quando o vento aumenta (vento espalha a poluição). Também espero diferença de ozônio entre os meses mais quentes e os mais amenos.

Código

dados <- airquality |>
  mutate(
    Mes = factor(Month,
                 levels = 5:9,
                 labels = c("Maio", "Junho", "Julho", "Agosto", "Setembro")),
    Faixa_Temp = cut(Temp,
                     breaks = c(-Inf, 70, 85, Inf),
                     labels = c("Fria", "Amena", "Quente"))
  )

glimpse(dados)

Rows: 153
Columns: 8
$ Ozone      <int> 41, 36, 12, 18, NA, 28, 23, 19, 8, NA, 7, 16, 11, 14, 18, 1…
$ Solar.R    <int> 190, 118, 149, 313, NA, NA, 299, 99, 19, 194, NA, 256, 290,…
$ Wind       <dbl> 7.4, 8.0, 12.6, 11.5, 14.3, 14.9, 8.6, 13.8, 20.1, 8.6, 6.9…
$ Temp       <int> 67, 72, 74, 62, 56, 66, 65, 59, 61, 69, 74, 69, 66, 68, 58,…
$ Month      <int> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,…
$ Day        <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
$ Mes        <fct> Maio, Maio, Maio, Maio, Maio, Maio, Maio, Maio, Maio, Maio,…
$ Faixa_Temp <fct> Fria, Amena, Amena, Fria, Fria, Fria, Fria, Fria, Fria, Fri…

2. Descrição estatística (gtsummary)

Código

dados |>
  select(Ozone, Solar.R, Wind, Temp, Mes, Faixa_Temp) |>
  tbl_summary(
    statistic = list(all_continuous() ~ "{median} ({p25}, {p75})"),
    missing_text = "Faltantes (NA)"
  )

Tabela 1

Characteristic	N = 153¹
Ozone	32 (18, 64)
Faltantes (NA)	37
Solar.R	205 (115, 259)
Faltantes (NA)	7
Wind	9.7 (7.4, 11.5)
Temp	79 (72, 85)
Mes
Maio	31 (20%)
Junho	30 (20%)
Julho	31 (20%)
Agosto	31 (20%)
Setembro	30 (20%)
Faixa_Temp
Fria	33 (22%)
Amena	86 (56%)
Quente	34 (22%)
¹ Median (Q1, Q3); n (%)

A tabela já mostra a mediana e os quartis das numéricas e a contagem das categóricas. Dá pra ver de cara que Ozone tem bastante NA — isso vai importar mais pra frente.

3. Gráfico de calor (heatmap de correlação)

Código

num <- dados |> select(Ozone, Solar.R, Wind, Temp)

# correlação só com as linhas completas (use = "complete.obs")
matriz_cor <- cor(num, use = "complete.obs")

matriz_long <- as.data.frame(as.table(matriz_cor))
colnames(matriz_long) <- c("var1", "var2", "correlacao")

ggplot(matriz_long, aes(var1, var2, fill = correlacao)) +
  geom_tile(color = "white") +
  geom_text(aes(label = round(correlacao, 2)), size = 4) +
  scale_fill_gradient2(low = "#2166ac", mid = "white", high = "#b2182b",
                       midpoint = 0, limits = c(-1, 1)) +
  labs(title = "Correlação entre as variáveis numéricas",
       x = NULL, y = NULL, fill = "Correlação") +
  theme_minimal()

Lendo o gráfico: as variáveis mais correlacionadas com o meu Ozone são a Temp (correlação positiva forte, perto de 0,7) e o Wind (correlação negativa, em torno de -0,6). Ou seja: dia mais quente, mais ozônio; dia mais ventoso, menos ozônio. O Solar.R tem uma correlação positiva fraca. Isso bate com o que eu esperava na seção 1.

4. Normalidade das variáveis

O que é uma distribuição normal

A distribuição normal é aquela curva em formato de sino, simétrica em volta da média. Nela, a média, a mediana e a moda ficam praticamente no mesmo ponto, e a maior parte dos dados se concentra perto do centro, ficando cada vez mais raro conforme se afasta. Muita coisa na estatística (vários testes) assume que os dados seguem mais ou menos esse formato, por isso vale a pena checar.

Histograma + densidade (3 variáveis)

Escolhi Ozone, Temp e Wind. Pro número de bins eu usei a regra da raiz quadrada do número de observações (sqrt(n)), que pra essa base dá em torno de 12 bins — é um número que mostra o formato sem ficar nem muito serrilhado nem muito grosseiro.

Código

n_bins <- round(sqrt(nrow(dados)))   # ~12

plot_hist <- function(df, coluna) {
  ggplot(df, aes(x = .data[[coluna]])) +
    geom_histogram(aes(y = after_stat(density)), bins = n_bins,
                   fill = "#80b1d3", color = "white") +
    geom_density(color = "#d6604d", linewidth = 1) +
    labs(title = paste("Histograma + densidade:", coluna), x = coluna, y = "Densidade") +
    theme_minimal()
}

plot_hist(dados, "Ozone")

Código

plot_hist(dados, "Temp")

Código

plot_hist(dados, "Wind")

Gráficos Q-Q

Código

plot_qq <- function(df, coluna) {
  ggplot(df, aes(sample = .data[[coluna]])) +
    stat_qq(color = "#4393c3") +
    stat_qq_line(color = "#d6604d") +
    labs(title = paste("Q-Q:", coluna), x = "Teórico", y = "Observado") +
    theme_minimal()
}

plot_qq(dados, "Ozone")

Código

plot_qq(dados, "Temp")

Código

plot_qq(dados, "Wind")

Teste de Shapiro-Wilk

Código

shapiro.test(dados$Ozone)


    Shapiro-Wilk normality test

data:  dados$Ozone
W = 0.87867, p-value = 2.79e-08

Código

shapiro.test(dados$Temp)


    Shapiro-Wilk normality test

data:  dados$Temp
W = 0.97617, p-value = 0.009319

Código

shapiro.test(dados$Wind)


    Shapiro-Wilk normality test

data:  dados$Wind
W = 0.98575, p-value = 0.1178

Conclusão sobre normalidade

Juntando tudo: o Ozone claramente não é normal — o histograma é bem puxado pra direita (assimétrico), os pontos do Q-Q fogem da linha nas pontas e o Shapiro dá p-valor bem abaixo de 0,05. O Wind está mais próximo de normal (Q-Q quase em cima da linha), e o Temp fica no meio do caminho. Pelo teste de Shapiro, a que mais se aproxima de uma distribuição normal é o Wind; Ozone é a mais longe de normal.

5. Cinco hipóteses

H1 — Quanto maior a temperatura, maior o ozônio (numérica x numérica)

Código

ggplot(dados, aes(Temp, Ozone)) +
  geom_point(color = "#3182bd") +
  geom_smooth(method = "lm", se = FALSE, color = "#d6604d") +
  labs(title = "Ozônio x Temperatura") + theme_minimal()

Código

cor.test(dados$Temp, dados$Ozone)


    Pearson's product-moment correlation

data:  dados$Temp and dados$Ozone
t = 10.418, df = 114, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5913340 0.7812111
sample estimates:
      cor 
0.6983603

O gráfico sobe da esquerda pra direita e o cor.test dá correlação positiva forte com p-valor bem pequeno. Confirma a hipótese: dia mais quente tende a ter mais ozônio.

H2 — O ozônio varia entre os meses (numérica x categórica)

Código

ggplot(dados, aes(Mes, Ozone, fill = Mes)) +
  geom_boxplot() +
  labs(title = "Ozônio por mês") + theme_minimal() + theme(legend.position = "none")

Código

dados |> select(Ozone, Mes) |>
  tbl_summary(by = Mes, statistic = all_continuous() ~ "{median} ({min}, {max})")

Characteristic	Maio N = 31¹	Junho N = 30¹	Julho N = 31¹	Agosto N = 31¹	Setembro N = 30¹
Ozone	18 (1, 115)	23 (12, 71)	60 (7, 135)	52 (9, 168)	23 (7, 96)
Unknown	5	21	5	5	1
¹ Median (Min, Max)

Código

summary(aov(Ozone ~ Mes, data = dados))

             Df Sum Sq Mean Sq F value   Pr(>F)    
Mes           4  29438    7359   8.536 4.83e-06 ***
Residuals   111  95705     862                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
37 observations deleted due to missingness

Os boxplots de julho e agosto ficam mais altos, e a ANOVA dá p-valor pequeno. Confirma: tem diferença de ozônio entre os meses, com os do verão mais altos.

H3 — Quanto mais vento, menos ozônio (numérica x numérica)

Código

ggplot(dados, aes(Wind, Ozone)) +
  geom_point(color = "#3182bd") +
  geom_smooth(method = "lm", se = FALSE, color = "#d6604d") +
  labs(title = "Ozônio x Vento") + theme_minimal()

Código

cor.test(dados$Wind, dados$Ozone)


    Pearson's product-moment correlation

data:  dados$Wind and dados$Ozone
t = -8.0401, df = 114, p-value = 9.272e-13
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.7063918 -0.4708713
sample estimates:
       cor 
-0.6015465

A nuvem de pontos desce e a correlação é negativa com p-valor pequeno. Confirma: mais vento, menos ozônio (o vento dispersa a poluição).

H4 — O ozônio é maior nas faixas de temperatura mais quentes (numérica x categórica)

Código

ggplot(dados, aes(Faixa_Temp, Ozone, fill = Faixa_Temp)) +
  geom_boxplot() +
  labs(title = "Ozônio por faixa de temperatura") +
  theme_minimal() + theme(legend.position = "none")

Código

summary(aov(Ozone ~ Faixa_Temp, data = dados))

             Df Sum Sq Mean Sq F value   Pr(>F)    
Faixa_Temp    2  54445   27222   43.51 9.73e-15 ***
Residuals   113  70698     626                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
37 observations deleted due to missingness

A faixa “Quente” tem a mediana de ozônio bem mais alta e a ANOVA confirma com p-valor pequeno. Confirma a hipótese.

H5 — Mais radiação solar, mais ozônio (numérica x numérica)

Código

ggplot(dados, aes(Solar.R, Ozone)) +
  geom_point(color = "#3182bd") +
  geom_smooth(method = "lm", se = FALSE, color = "#d6604d") +
  labs(title = "Ozônio x Radiação Solar") + theme_minimal()

Código

cor.test(dados$Solar.R, dados$Ozone)


    Pearson's product-moment correlation

data:  dados$Solar.R and dados$Ozone
t = 3.8798, df = 109, p-value = 0.0001793
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.173194 0.502132
sample estimates:
      cor 
0.3483417

Aqui a relação existe mas é mais fraca: a correlação é positiva e o p-valor é significativo, só que o efeito é bem menor que o da temperatura. Confirma parcialmente: influencia, mas pouco.

6. Completude de cada variável

Código

dados |>
  summarise(across(everything(), ~ mean(!is.na(.x)) * 100)) |>
  pivot_longer(everything(), names_to = "variavel", values_to = "completude_pct") |>
  arrange(completude_pct)

# A tibble: 8 × 2
  variavel   completude_pct
  <chr>               <dbl>
1 Ozone                75.8
2 Solar.R              95.4
3 Wind                100  
4 Temp                100  
5 Month               100  
6 Day                 100  
7 Mes                 100  
8 Faixa_Temp          100

A maioria das variáveis está 100% completa. As que têm buraco são o Ozone (a com mais NA) e o Solar.R (alguns NA).

7. Imputação dos dados faltantes

Vou imputar Ozone e Solar.R pela mediana. Escolhi a mediana (e não a média) porque essas duas variáveis são assimétricas e têm valores altos puxando a média pra cima — a mediana é mais robusta a esses valores extremos, então não distorce tanto os dados. É uma imputação simples, que mantém o tamanho da base sem inventar valores muito fora da realidade.

Código

dados_imp <- dados |>
  mutate(
    Ozone   = ifelse(is.na(Ozone),   median(Ozone,   na.rm = TRUE), Ozone),
    Solar.R = ifelse(is.na(Solar.R), median(Solar.R, na.rm = TRUE), Solar.R)
  )

# conferindo que não sobrou NA
colSums(is.na(dados_imp))

     Ozone    Solar.R       Wind       Temp      Month        Day        Mes 
         0          0          0          0          0          0          0 
Faixa_Temp 
         0

8. Modelos de regressão linear

Modelo 1: ozônio explicado só pela temperatura. Modelo 2: a mesma coisa, mais o vento.

Código

modelo1 <- lm(Ozone ~ Temp, data = dados_imp)
modelo2 <- lm(Ozone ~ Temp + Wind, data = dados_imp)

summary(modelo1)


Call:
lm(formula = Ozone ~ Temp, data = dados_imp)

Residuals:
    Min      1Q  Median      3Q     Max 
-36.464 -16.399  -3.866  10.314 122.691 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -104.0802    15.6655  -6.644 5.21e-10 ***
Temp           1.8443     0.1997   9.236  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 23.3 on 151 degrees of freedom
Multiple R-squared:  0.361, Adjusted R-squared:  0.3568 
F-statistic: 85.31 on 1 and 151 DF,  p-value: < 2.2e-16

Código

summary(modelo2)


Call:
lm(formula = Ozone ~ Temp + Wind, data = dados_imp)

Residuals:
    Min      1Q  Median      3Q     Max 
-36.125 -13.914  -3.643  12.889 106.546 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -41.8319    19.6694  -2.127   0.0351 *  
Temp          1.3876     0.2102   6.603 6.56e-10 ***
Wind         -2.6792     0.5646  -4.745 4.81e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.8 on 150 degrees of freedom
Multiple R-squared:  0.4444,    Adjusted R-squared:  0.437 
F-statistic: 59.99 on 2 and 150 DF,  p-value: < 2.2e-16

Código

# comparando pelo R² ajustado e pelo AIC
tibble(
  modelo = c("modelo1 (Temp)", "modelo2 (Temp + Wind)"),
  r2_ajustado = c(summary(modelo1)$adj.r.squared, summary(modelo2)$adj.r.squared),
  AIC = c(AIC(modelo1), AIC(modelo2))
)

# A tibble: 2 × 3
  modelo                r2_ajustado   AIC
  <chr>                       <dbl> <dbl>
1 modelo1 (Temp)              0.357 1402.
2 modelo2 (Temp + Wind)       0.437 1382.

Escolha: fico com o modelo 2. O R² ajustado dele é maior (explica mais a variação do ozônio) e o AIC é menor (penaliza menos o ajuste), então adicionar o vento valeu a pena — ele traz informação nova que a temperatura sozinha não tinha.

Testando o modelo final

Código

# previsão pra um dia quente e com pouco vento
novo_dia <- data.frame(Temp = 85, Wind = 5)
predict(modelo2, newdata = novo_dia)

       1 
62.71754

Código

# resíduos do modelo escolhido
ggplot(data.frame(ajustado = fitted(modelo2), residuo = resid(modelo2)),
       aes(ajustado, residuo)) +
  geom_point(color = "#3182bd") +
  geom_hline(yintercept = 0, color = "#d6604d", linetype = "dashed") +
  labs(title = "Resíduos x valores ajustados (modelo 2)",
       x = "Ajustado", y = "Resíduo") + theme_minimal()

A previsão pra um dia quente (85°F) e pouco vento dá um ozônio alto, que é exatamente o que as hipóteses anteriores apontaram. O gráfico de resíduos não mostra um padrão muito grave, então o modelo está aceitável pro nível desse trabalho.

9. Fontes e referências

Documentação da base airquality: https://stat.ethz.ch/R-manual/R-devel/library/datasets/html/airquality.html
Pacote gtsummary: https://www.danieldsjoberg.com/gtsummary/
Pacote ggplot2: https://ggplot2.tidyverse.org/
Quarto: https://quarto.org/