Pollution de l’air

La pollution de l’air est un problème de santé publique majeur. De nombreuses études ont démontré l’influence sur la santé de certains composés chimiques, dont l’ozone (\(O_{3}\)), tout particulièrement sur les personnes sensibles (nouveaux-nés, asthmatiques, personnes âgées). Il est donc important de savoir prévoir les pics de concentration de l’ozone. On sait que la concentration en ozone varie avec la température : plus la température est élevée, plus la concentration en ozone est importante. Cette relation vague doit être précisée en vue de prévoir les pics d’ozone. Dans ce but, l’association Air Breizh (surveillance de la qualité de l’air en Bretagne) mesure depuis 1994 la concentration en \(O_3\) (en \(\mu g/L\)) toutes les dix minutes et obtient donc le maximum journalier de la concentration en \(O_{3}\) , que l’on note désormais \(O3\). Air Breizh collecte également des données météorologiques correspondant à ces mesures d’ozone, dont la température à 12h, que l’on note \(T12\). À partir de ces données (rassemblées dans le fichier data/ozone_simple.txt), on cherche à expliquer la concentration maximale \(03\) à l’aide de la température T12.

Préliminaire : analyse descriptive

Importer les données ozone_simple.txt.

library(tidyverse)
library(plotly)
data <- readr::read_delim(
  here::here("data", "ozone_simple.txt"), delim = ";", 
  col_types = c(col_double(), col_double())
)

Calculer la moyenne, la médiane, la variance de la concentration maximale 03 et de la température T12.

data %>%
  summarise(
    mean03 = mean(O3),
    meanT12 = mean(T12),
    medianO3 = median(O3),
    medianT12 = median(T12),
    var03 = var(O3),
    varT12 = var(T12),
    maxO3 = max(O3),
    maxT12 = max(T12)
  ) %>%
    t()

##                [,1]
## mean03     86.30000
## meanT12    20.32000
## medianO3   83.90000
## medianT12  19.75000
## var03     571.24939
## varT12     21.85224
## maxO3     139.00000
## maxT12     29.50000

Visualisation

Tracer l’histogramme de la concentration maximale 03, puis celui de la température T12.

source(here::here("utils", "plotting-funcs.R"))

## here() starts at /home/gml/Documents/Master/M1_BIBS/S2/MultiVar/TPs

p.density.o3 <- data %>%
    ggplot(aes(x = O3)) +
      geom_histogram(
        aes(y = stat(density)), colour="black", fill="white", 
        binwidth = optim.binwidth(data$O3)
      ) + 
        geom_density(alpha=.2, fill="blue")


p.density.t12 <- data %>%
  ggplot(aes(x = T12)) +
    geom_histogram(
        aes(y = stat(density)), colour="black", fill="white", 
        binwidth = optim.binwidth(data$T12)
      ) + 
        geom_density(alpha=.2, fill="blue")

ggplotly(p.density.o3)

ggplotly(p.density.t12)

Tracer, sous la forme d’un nuage de point, le graphe de la concentration maximale 03 en fonction de la température T12.

plot.scatter <- data %>%
  ggplot(aes(x = T12, y = O3)) +
    geom_point()

ggplotly(plot.scatter)

Calcul des estimateurs des moindres carrés

Écrire un modèle de régression linéaire permettant d’expliquer la concentration maximale 03 en fonction de la température T12.

\[ O3 = \beta_{1} + \beta_{2}T12 + \epsilon \]

Utiliser R pour estimer les paramètres de ce modèle :

D’abord en utilisant les formules du cours,
ensuite à l’aide de la fonction lm.
Donner les estimations de l’ordonnée à l’origine, la pente et la variance résiduelle.

beta2_chap <- cov(data$T12, data$O3) / var(data$T12)
beta1_chap <- mean(data$O3) - beta2_chap*mean(data$T12)
glue::glue("Estimation manuelle\n beta1 ~ {beta1_chap}, beta2 ~ {beta2_chap}")

## Estimation manuelle
## beta1 ~ 31.4149732899996, beta2 ~ 2.70103477903545

modele <- lm(O3 ~ T12, data = data)
beta1_lm <- modele$coefficients[["(Intercept)"]]
beta2_lm <- modele$coefficients[["T12"]]
var.resid <- var(modele$residuals)
glue::glue("Estimation avec lm\n beta1 ~ {beta1_lm}, beta2 ~ {beta2_lm}, var.res ~ {var.resid}")

## Estimation avec lm
## beta1 ~ 31.4149732899996, beta2 ~ 2.70103477903545, var.res ~ 411.824392927854

Visualisation de la droite de régression

Tracer la droite de régression (superposée avec le nuage de points représentant les observations).

plot.regline <- data %>%
  ggplot(mapping = aes(x=T12, y=O3)) + 
    geom_point() +
      geom_abline(slope = beta2_lm, intercept = beta1_lm, colour="blue")
ggplotly(plot.regline)

Analyse des résidus

Faire un graphe des résidus. Que constate-t-on ?

residus <- tibble(T12=data$T12, residuals=modele$residuals)

my.qqplot <- residus %>%
  ggplot(mapping = aes(sample = residuals)) +
    stat_qq() + 
      stat_qq_line()
ggplotly(my.qqplot)

Les résidus ne semblent d’être i.i.d. Leur distribution semble d’être dépendante de la valeur de la température \(T12\) ce qui pourrait invalider le modèle qui est déjà mauvais en termes de qualité prédictive \(r^2 = 0.2790\).

L’estimation de la densité des résidus à un léger caractère bimodale, ce qui pourrait violer l’hypothèse de normalité des résidus. Or, le shapiro.test montre que l’échantillon ne supporte suffisamment le rejet de l’hypothèse de nulle de normalité à un \(95\%\) de confiance.

plot.reg.res <- residus %>%
  ggplot(mapping = aes(x=T12, y=residuals)) +
    geom_point() + 
      geom_smooth()

plot.reg.res.dist <- residus %>%
  ggplot(aes(x = residuals)) +
    geom_histogram(
        aes(y = stat(density)), colour="black", fill="white", 
        binwidth = optim.binwidth(residus$residuals)
      ) + 
        geom_density(alpha=.2, fill="blue")

ggplotly(plot.reg.res)

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

ggplotly(plot.reg.res.dist)

shapiro.test(modele$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modele$residuals
## W = 0.97852, p-value = 0.4912

TP1 Régression linéaire simple : pollution de l’air

Gustavo Magaña López