Commit 089e687d authored by Laurence Viry's avatar Laurence Viry
Browse files

correct ozone

parent aa6edcb4
Pipeline #11227 passed with stages
in 1 minute and 20 seconds
This diff is collapsed.
%% Cell type:markdown id: tags:
## Etude de la polution de l'air en fonction de données climatiques
(*Cet exemple est tiré du livre "Statistiques avec R" de Pierre-André Cornillon (Presse Universitaire de Rennes), il peut être intéressant pour ceux qui sont intéressés par l'interprétation des résultats de le consulter*).
La pollution de l'air est actuellement une des préoccupations majeures de santé publique. Des associations de surveillance de la qualité de l'air existent sur tout le territoire français et mesurent la **concentration des polluants** ainsi que **les conditions météorologiques** comme la température, la nébulosité, le vent,$\ldots$
......@@ -22,14 +19,14 @@
%% Cell type:code id: tags:
``` R
# Lecture des données
setwd('/home/viryl/notebooks/CED-IntroR/TP/data')
...
ozone <- read.csv("ozone.txt",sep="",row.names=1,header=TRUE)
names(ozone)
# Extraction des données utiles
ozone1 <-
ozone1 <- ozone[,c("maxO3","T12")]
names(ozone1)
```
%% Cell type:markdown id: tags:
......@@ -37,11 +34,11 @@
%% Cell type:code id: tags:
``` R
# Statistiques descriptives des données de l'étude
...
summary(ozone1)
```
%% Cell type:markdown id: tags:
#### Analyse du lien entre ces deux variables par une représentation graphique.
......@@ -50,11 +47,11 @@
%% Cell type:code id: tags:
``` R
# Nuage de point ("maxO3","T12")
...
plot(maxO3~T12,data=ozone1,pch=15,cex=0.5)
```
%% Cell type:markdown id: tags:
#### Estimation des paramètres du modèle de régression simple.<br\>
......@@ -72,15 +69,15 @@
```
%% Cell type:code id: tags:
``` R
reg.simple <-
reg.simple <- lm(maxO3~T12,data=ozone1)
names(reg.simple)
reg.simple$coefficients
class(reg.simple)
...
summary(reg.simple)
```
%% Cell type:markdown id: tags:
#### Tracer la droite de régression sur le nuage de point obtenue en 3.
......@@ -90,20 +87,21 @@
%% Cell type:code id: tags:
``` R
# Utilisation de la fonction abline
...
plot(maxO3~T12,data=ozone1,pch=15,cex=0.5)
abline(reg.simple,col=2)
```
%% Cell type:code id: tags:
``` R
plot(maxO3~T12,data=ozone1,pch=15,cex=0.5)
grillex <- ...
grilley <- ...
...
grillex <- seq(min(ozone1[,"T12"]),max(ozone1[,"T12"]),length=100)
grilley <- reg.simple$coef[1]+ reg.simple$coef[2]*grillex
lines(grillex,grilley,col=2)
```
%% Cell type:markdown id: tags:
#### Analyser les résidus
......@@ -120,11 +118,13 @@
%% Cell type:code id: tags:
``` R
xnew <- 19
...
xnew <- as.data.frame(xnew)
colnames(xnew) <- "T12"
predict(reg.simple,xnew, interval="pred")
```
%% Cell type:markdown id: tags:
Pour représenter sur un même graphe l'intervalle de confiance d'une valeur lissée et l'intervalle de confiance d'une valeur prédite, nous calculons ces intervales sur l'ensemble des points ayant servi à la construction de la droite de régression.
......@@ -144,11 +144,11 @@
%% Cell type:markdown id: tags:
### Regression multiple
La régression linéaire multiple permet d'expliquer et/ou prédire une variable quantitative par p variables quantitatives, c'est une généralisation du modèle de régression simple.
La régression linéaire multiple permet d'expliquer et/ou prédire une variable quantitative par p variables qunatitatives, c'est une généralisation du modèle de régression simple.
On revient à notre problème et on veut analyser la relation entre le maximum journalier de la concentration en ozone et la température à différentes heures de la journée, la nébulosité à différentes heures de la journée, la projection du vent sur l'axe Est-Ouest à différentes heures de la journée et la concentration maximale de la veille du jour considéré.
** Les diférentes étapes **
......@@ -158,11 +158,11 @@
``` R
# Lecture des données
names(ozone)
# Extraction des données utiles
ozone2 <- ...
ozone2 <- ozone[,c('maxO3','T9','T12','T15','Ne9','Ne12','Ne15','Vx9','Vx12','Vx15','maxO3v')]
names(ozone2)
summary(ozone2)
```
%% Cell type:markdown id: tags:
......@@ -171,30 +171,21 @@
Pour valider les données, on effectue une analyse univariée des variables(statistiques descriptives, histogramme,...)
Le nombre de variables n'étant pas très élevé, on peut représenter les variables deux par deux sur un même graphique (fonction **pairs**). On pourrait également explorer les données avec une ACP (package FactoMineR, FactoShiny).
%% Cell type:code id: tags:
``` R
# plot
```
%% Cell type:markdown id: tags:
#### Estimer les paramètres
Pour estimer les paramètres, il faut écrire le modèle. On utilise la fonction **lm** avec une formule décrivant le modèle.<br\>
**maxO3 ~ T9+ T12+T15+Ne9+Ne12+Ne15+Vx9+Vx12+Vx15+maxO3v'** <br\>
%% Cell type:code id: tags:
``` R
# regression multiple
reg.mul <-
reg.mul <- lm(maxO3 ~ T9+ T12+T15+Ne9+Ne12+Ne15+Vx9+Vx12+Vx15+maxO3v, data=ozone2)
```
%% Cell type:code id: tags:
``` R
......@@ -207,18 +198,18 @@
#### Choix des variables
Il est possible avec l'analyse de **reg.mul** de faire un choix des variables à la main. On enlève la moins significative puis on recalcule les estimations et ainsi de suite.
Il existe aussi un package R qui traite du choix des variables: le package **leaps** et la fonction **regsubsets**
Il existe aussi un package R qui traite du choix des variables: le package **leaps**
%% Cell type:code id: tags:
``` R
library("leaps", lib.loc="~/R/lib")
library("lattice")
choix <-
choix <- regsubsets(maxO3 ~ ., data=ozone2, nbest=1,nvmax=11)
plot(choix)
```
%% Cell type:markdown id: tags:
......@@ -226,30 +217,28 @@
%% Cell type:code id: tags:
``` R
# regression après selection des variables
reg.fin <-
reg.fin <- lm(maxO3~T12+Ne9+Vx9+maxO3v,data=ozone2)
summary(reg.fin)
```
%% Cell type:markdown id: tags:
#### Analyser les résidus:
Même traitement qu'avec la régression simple.
#### Analyser les résidus: même traitement qu'avec la régression simple
#### Prévoir une nouvelle valeur: même traitement qu'avec la régression simple
%% Cell type:code id: tags:
``` R
# Prediction + interval de confiance
xnew <- matrix(c(19,3,30,6.5),nrow=1)
colnames(xnew) <- c("T12","Ne9","Vx9","maxO3v")
xnew <-
ynew <-
xnew <- as.data.frame(xnew)
ynew <- predict(reg.fin,xnew,interval="pred")
class(ynew)
dimnames(ynew)
```
%% Cell type:markdown id: tags:
......@@ -268,11 +257,11 @@
``` R
# Lecture des données
names(ozone)
# Extraction des données utiles
ozone3 <- ...
ozone3 <- ozone[,c('maxO3','vent')]
names(ozone3)
summary(ozone3)
```
%% Cell type:markdown id: tags:
......@@ -285,11 +274,11 @@
%% Cell type:code id: tags:
``` R
# boxplot par modalité
...
plot(maxO3 ~vent, data=ozone3,pch=15,cex=0.5)
```
%% Cell type:markdown id: tags:
#### Analyse de la significativité du facteur
......@@ -298,18 +287,18 @@
%% Cell type:code id: tags:
``` R
# Estimation des parametres
reg.aov1 <-
reg.aov1 <- lm(maxO3 ~vent,data=ozone3)
```
%% Cell type:code id: tags:
``` R
# Tableau d'analyse de variance
anova(reg.aov1)
```
%% Cell type:markdown id: tags:
#### Analyser les résidus
......@@ -359,11 +348,11 @@
``` R
# Lecture des données
names(ozone)
# Extraction des données utiles
ozone4 <- ...
ozone4 <- ozone[,c('maxO3','vent','pluie')]
names(ozone4)
summary(ozone4)
```
%% Cell type:markdown id: tags:
......@@ -376,11 +365,11 @@
%% Cell type:code id: tags:
``` R
# Representation des donnees
...
boxplot(maxO3~vent*pluie,data=ozone4,cex=0.5)
# Interaction
par(mfrow=c(1,2))
with(ozone,interaction.plot(vent,pluie,maxO3))
with(ozone,interaction.plot(pluie,vent,maxO3))
```
......@@ -392,19 +381,19 @@
%% Cell type:code id: tags:
``` R
# Choisir le modèle
mod.int <-
mod.int <- lm(maxO3~vent*pluie,data=ozone4)
anova(mod.int)
```
%% Cell type:code id: tags:
``` R
# Choisir le modèle
mod.ssint <-
mod.ssint <- lm(maxO3~vent+pluie,data=ozone4)
anova(mod.ssint)
```
%% Cell type:markdown id: tags:
......@@ -413,11 +402,11 @@
Comme pour les méthodes précédentes, on utilisera la fonction **summary** pour aider l'interprétation.
%% Cell type:code id: tags:
``` R
...
summary(mod.ssint)
```
%% Cell type:markdown id: tags:
Une interprétation plus approfondie serait nécessaire mais ça n'est pas le propos de ce cours.
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment