Prévision météorologique sur la Ville de Tours
Charles Vitry
Clovis Deletre
Régression Linéaire Multiple
Présentation des données
Nous étions à la recherche d’un jeu de données de minimum 300 observations et 15 variables quantitatives dont une à expliquer, Nous avions des jeux de données provenant de sondages, mais leurs données étaient qualitative, donc non exploitable dans notre contexte d’étude.
Après recherche, le jeu de données choisis est météorologique, nous souhaitons prévoir la température du lendemain en fonction des autres variables quantitative. Ce jeu de données a été obtenu dans les open data de la ville de Tours, le fichier est trouvable à l’adresse suivante : https://data.tours-metropole.fr/explore/dataset/observation-meteorologique-historiques-tours-synop/export/?sort=date.
# Import de la base de données
On ajoute le paramètre UTF 8 car les noms de variables sont en français.
#Windows only
# data <-
# read.csv("~/PrevisionMeteoTours/observation-meteorologique-historiques-tours-synop.csv",
# encoding="UTF-8",
# sep=";")
#Linux & Windows
data <-
read.csv("observation-meteorologique-historiques-tours-synop.csv",
encoding="UTF-8",
sep=";")D’après l’hébergeur les paramètres atmosphériques sont mesurés (température, humidité, direction et force du vent, pression atmosphérique…) ou observés (description des nuages, visibilité…) depuis la surface terrestre.
Après lecture de la colonne Date, les données sont collectés de 2010 jusqu’à maintenant, il en résulte 34854 observations.
Nous avons cependant une limite à notre analyse. Puisque l’ensemble de ces observations est tirés d’une seule station d’observation. 
Nous supprimons alors les collonnes correspondant à des indications géographiques
#numero de la station météo
data$ID.OMM.station <- NULL
data$ID.OMM.station<- NULL
#coords géographique
data$Coordonnees <- NULL
data$Coordonnees <- NULL
data$Latitude <- NULL
data$Latitude <- NULL
data$Longitude <- NULL
data$Longitude <- NULL
#Nom de la ville et département
data$Nom <- NULL
data$Nom <- NULL
data$communes..code. <- NULL
data$communes..code. <- NULL
data$communes..name. <- NULL
data$communes..name. <- NULL
data$EPCI..name. <- NULL
data$EPCI..name. <- NULL
data$EPCI..code. <- NULL
data$EPCI..code. <- NULL
data$department..name. <- NULL
data$department..name. <- NULL
data$department..code. <- NULL
data$department..code. <- NULL
data$region..name. <- NULL
data$region..name. <- NULL
data$region..code.<- NULL
data$region..code. <- NULL
Les colonnes vides ou constantes sont inutilisables
data$Niveau.baromÃ.trique <- NULL
data$Niveau.barométrique <- NULL
data$TempÃ.rature.minimale.sur.24.heures <- NULL
data$Température.minimale.sur.24.heures...C. <- NULL
data$TempÃ.rature.minimale.sur.24.heures..Â.C. <- NULL
data$MÃ.thode.de.mesure.TempÃ.rature.du.thermomÃ.tre.mouillÃ. <- NULL
data$Méthode.de.mesure.Température.du.thermomètre.mouillé <- NULL
data$TempÃ.rature.du.thermomÃ.tre.mouillÃ. <- NULL
data$Température.du.thermomètre.mouillé <- NULL
data$Hauteur.totale.de.la.couche.de.neige..glace..autre.au.sol <- NULL
data$Hauteur.totale.de.la.couche.de.neige..glace..autre.au.sol <- NULL
data$NÃ.bulositÃ..couche.nuageuse.4 <- NULL
data$Nébulosité.couche.nuageuse.4 <- NULL
data$Hauteur.de.base.4 <- NULL
data$Hauteur.de.base.4 <- NULL
data$TempÃ.rature.minimale.sur.24.heures <- NULL
data$Température.minimale.sur.24.heures <- NULL
data$TempÃ.rature.maximale.sur.24.heures <- NULL
data$Température.maximale.sur.24.heures <- NULL
data$Altitude <- NULL
data$Altitude <- NULL
data$NÃ.bulositÃ..couche.nuageuse.1 <- NULL
data$Nébulosité.couche.nuageuse.1 <- NULL
data$Etat.du.sol <- NULL
data$Etat.du.sol <- NULL
data$Periode.de.mesure.de.la.rafale <- NULL
data$Periode.de.mesure.de.la.rafale <- NULL
data$NÃ.bulositÃ...des.nuages.de.l..Ã.tage.infÃ.rieur <- NULL
data$Nébulosité..des.nuages.de.l..étage.inférieur <- NULL
data$TempÃ.rature.minimale.du.sol.sur.12.heures..en.Â.C. <- NULL
data$Température.minimale.du.sol.sur.12.heures..en..C. <- NULL
Retirons les variables qualitatives.
Type de tendance barométrique et temps présent est mis quali- mi quanti, nous ne les traiterons pas
data$Type.de.tendance.baromÃ.trique <- NULL
data$Type.de.tendance.barométrique <- NULL
data$Temps.prÃ.sent <- NULL
data$Temps.présent <- NULL
data$Temps.prÃ.sent.1 <- NULL
data$Temps.présent.1 <- NULL
data$Temps.passÃ..1 <- NULL
data$Temps.passé.1 <- NULL
data$Temps.passÃ..1.1 <- NULL
data$Temps.passé.1.1 <- NULL
data$Temps.passÃ..2 <- NULL
data$Temps.passé.2 <- NULL
data$Type.de.tendance.baromÃ.trique.1 <- NULL
data$Type.de.tendance.barométrique.1 <- NULL
data$mois_de_l_annee <- NULL
data$mois_de_l_annee <- NULL
Nous n’avons pas supprimé la variable date car elle nous est utile à trier chronologiquement les données.
data <- data[order(as.Date(data$Date)),]
Observons les valeurs manquantes des variables restantes
#1
#sum(is.na(data$Pression.au.niveau.mer))
#2
#sum(is.na(data$Variation.de.pression.en.3.heures))
#3
#sum(is.na(data$Type.de.tendance.baromÃ.trique))
#4
#sum(is.na(data$Direction.du.vent.moyen.10.mn))
#5
#sum(is.na(data$Vitesse.du.vent.moyen.10.mn))
#6
#sum(is.na(data$TempÃ.rature))
#7
#sum(is.na(data$Point.de.rosÃ.e))
#8
#sum(is.na(data$HumiditÃ.))
#9
#sum(is.na(data$VisibilitÃ..horizontale))
#10
#sum(is.na(data$NebulositÃ..totale))
#11
#sum(is.na(data$NÃ.bulositÃ...des.nuages.de.l..Ã.tage.infÃ.rieur))
#12
#sum(is.na(data$Hauteur.de.la.base.des.nuages.de.l.Ã.tage.infÃ.rieur))
#13
#sum(is.na(data$Type.des.nuages.de.l.Ã.tage.infÃ.rieur))
#14
#sum(is.na(data$Type.des.nuages.de.l.Ã.tage.moyen))
#15
#sum(is.na(data$Type.des.nuages.de.l.Ã.tage.supÃ.rieur))
#16
#sum(is.na(data$Pression.station))
#17
#sum(is.na(data$GÃ.opotentiel))
#18
#sum(is.na(data$Variation.de.pression.en.24.heures))
#19
#sum(is.na(data$TempÃ.rature.minimale.sur.12.heures))
#20
#sum(is.na(data$TempÃ.rature.maximale.sur.12.heures))
#21
#sum(is.na(data$TempÃ.rature.minimale.du.sol.sur.12.heures))
#22
#sum(is.na(data$Rafale.sur.les.10.derniÃ.res.minutes))
#23
#sum(is.na(data$Rafales.sur.une.pÃ.riode))
#24
#sum(is.na(data$Periode.de.mesure.de.la.rafale))
#25
#sum(is.na(data$Etat.du.sol))
#26
#sum(is.na(data$Hauteur.de.la.neige.fraÃ.che))
#27
#sum(is.na(data$Periode.de.mesure.de.la.neige.fraiche))
#28
#sum(is.na(data$PrÃ.cipitations.dans.la.derniÃ.re.heure))
#29
#sum(is.na(data$PrÃ.cipitations.dans.les.3.derniÃ.res.heures))
#30
#sum(is.na(data$PrÃ.cipitations.dans.les.6.derniÃ.res.heures))
#31
#sum(is.na(data$PrÃ.cipitations.dans.les.12.derniÃ.res.heures))
#32
#sum(is.na(data$PrÃ.cipitations.dans.les.24.derniÃ.res.heures))
#33
#sum(is.na(data$PhÃ.nomÃ.ne.spÃ.cial.1))
#34
#sum(is.na(data$PhÃ.nomÃ.ne.spÃ.cial.2))
#35
#sum(is.na(data$PhÃ.nomÃ.ne.spÃ.cial.3))
#36
#sum(is.na(data$PhÃ.nomÃ.ne.spÃ.cial.4))
#37
#sum(is.na(data$NÃ.bulositÃ..couche.nuageuse.1))
#38
#sum(is.na(data$Type.nuage.1))
#39
#sum(is.na(data$Hauteur.de.base.1))
#40
#sum(is.na(data$NÃ.bulositÃ..couche.nuageuse.2))
#41
#sum(is.na(data$Type.nuage.2))
#42
#sum(is.na(data$NÃ.bulositÃ..couche.nuageuse.3))
#43
#sum(is.na(data$Type.nuage.3))
#44
#sum(is.na(data$Hauteur.de.base.3))
#45
#sum(is.na(data$Type.nuage.4))
#46
#sum(is.na(data$Temps.prÃ.sent.1))
#47
#sum(is.na(data$TempÃ.rature..Â.C.))
#48
#sum(is.na(data$TempÃ.rature.minimale.sur.12.heures..Â.C.))
#49
#sum(is.na(data$TempÃ.rature.maximale.sur.12.heures..Â.C.))
#50
#sum(is.na(data$TempÃ.rature.maximale.sur.24.heures..Â.C.))
#51
#sum(is.na(data$TempÃ.rature.minimale.du.sol.sur.12.heures..en.Â.C.))
#52
#sum(is.na(data$Altitude))
#53
#sum(is.na(data$mois_de_l_annee))On en déduit cela
plus 20 000 manquants:- $Type.des.nuages.de.l.Ã.tage.infÃ.rieur
- $Type.des.nuages.de.l.Ã.tage.moyen
- $Type.des.nuages.de.l.Ã.tage.supÃ.rieur
- $GÃ.opotentiel
- $TempÃ.rature.minimale.sur.12.heures
- $TempÃ.rature.maximale.sur.12.heures
- $Hauteur.de.la.neige.fraÃ.che
- $Periode.de.mesure.de.la.neige.fraiche
- $PhÃ.nomÃ.ne.spÃ.cial.1
- $PhÃ.nomÃ.ne.spÃ.cial.2
- $PhÃ.nomÃ.ne.spÃ.cial.3
- $PhÃ.nomÃ.ne.spÃ.cial.4
- $Type.nuage.1
- $Hauteur.de.base.1
- $NÃ.bulositÃ..couche.nuageuse.2
- $Type.nuage.2
- $NÃ.bulositÃ..couche.nuageuse.3
- $Type.nuage.3
- $Hauteur.de.base.3
- $Type.nuage.4
- $TempÃ.rature.minimale.sur.12.heures..Â.C.
- $TempÃ.rature.maximale.sur.12.heures..Â.C.
- $TempÃ.rature.maximale.sur.24.heures..Â.C.
Plus de 10 000:
- $NebulositÃ..totale
- $Hauteur.de.la.base.des.nuages.de.l.Ã.tage.infÃ.rieur
- $Variation.de.pression.en.24.heures
- $TempÃ.rature.minimale.du.sol.sur.12.heures
- $Rafale.sur.les.10.derniÃ.res.minutes
- $Etat.du.sol
- $PrÃ.cipitations.dans.les.6.derniÃ.res.heures
- $PrÃ.cipitations.dans.les.12.derniÃ.res.heures
- $PrÃ.cipitations.dans.les.24.derniÃ.res.heures
- $NÃ.bulositÃ..couche.nuageuse.1
- $TempÃ.rature.minimale.du.sol.sur.12.heures..en.Â.C.
- $NÃ.bulositÃ…des.nuages.de.l..Ã.tage.infÃ.rieur
Création d’un jeu de donnée où on supprimme que les plus de 30 000
#data2 <- data
#data2$GÃ.opotentiel <- NULL
#data2$TempÃ.rature.maximale.sur.24.heures..Â.C. <- NULL
#data2$Type.nuage.4 <- NULL
#data2$PhÃ.nomÃ.ne.spÃ.cial.4 <- NULL
#data2$Type.nuage.2 <- NULL
#data2$NÃ.bulositÃ..couche.nuageuse.3 <- NULL
#data2$Type.nuage.3 <- NULL
#data2$Hauteur.de.base.3 <- NULL
#data2$Type.nuage.4 <- NULLTest des régressions sur cet ensemble : exactement les mêmes résultats.
Suppression des plus de 20 000
data$Type.des.nuages.de.l.Ã.tage.infÃ.rieur <- NULL
data$Type.des.nuages.de.l.étage.inférieur <- NULL
data$Type.des.nuages.de.l.Ã.tage.moyen <- NULL
data$Type.des.nuages.de.l.étage.moyen <- NULL
data$Type.des.nuages.de.l.Ã.tage.supÃ.rieur <- NULL
data$Type.des.nuages.de.l.étage.supérieur <- NULL
data$GÃ.opotentiel <- NULL
data$Géopotentiel <- NULL
data$TempÃ.rature.minimale.sur.12.heures <- NULL
data$Température.minimale.sur.12.heures <- NULL
data$TempÃ.rature.maximale.sur.12.heures <- NULL
data$Température.maximale.sur.12.heures <- NULL
data$Hauteur.de.la.neige.fraÃ.che <- NULL
data$Hauteur.de.la.neige.fraîche <- NULL
data$Periode.de.mesure.de.la.neige.fraiche <- NULL
data$Periode.de.mesure.de.la.neige.fraiche <- NULL
data$PhÃ.nomÃ.ne.spÃ.cial.1 <- NULL
data$Phénomène.spécial.1 <- NULL
data$PhÃ.nomÃ.ne.spÃ.cial.2 <- NULL
data$Phénomène.spécial.2 <- NULL
data$PhÃ.nomÃ.ne.spÃ.cial.3 <- NULL
data$Phénomène.spécial.3 <- NULL
data$PhÃ.nomÃ.ne.spÃ.cial.4 <- NULL
data$Phénomène.spécial.4 <- NULL
data$Type.nuage.1 <- NULL
data$Type.nuage.1 <- NULL
data$Hauteur.de.base.1 <- NULL
data$Hauteur.de.base.1 <- NULL
data$NÃ.bulositÃ..couche.nuageuse.2 <- NULL
data$Nébulosité.couche.nuageuse.2 <- NULL
data$Type.nuage.2 <- NULL
data$Type.nuage.2 <- NULL
data$NÃ.bulositÃ..couche.nuageuse.3 <- NULL
data$Nébulosité.couche.nuageuse.3 <- NULL
data$Type.nuage.3 <- NULL
data$Type.nuage.3 <- NULL
data$Hauteur.de.base.3 <- NULL
data$Hauteur.de.base.3 <- NULL
data$Type.nuage.4 <- NULL
data$Type.nuage.4 <- NULL
data$TempÃ.rature.minimale.sur.12.heures..Â.C. <- NULL
data$Température.minimale.sur.12.heures...C. <- NULL
data$TempÃ.rature.maximale.sur.12.heures..Â.C. <- NULL
data$Température.maximale.sur.12.heures...C. <- NULL
data$TempÃ.rature.maximale.sur.24.heures..Â.C. <- NULL
data$Température.maximale.sur.24.heures...C. <- NULL
Il s’agit d’une conversion de Kelvin en celsius, cette ligne n’a pas d’utilité
data$TempÃ.rature..Â.C. <- NULL
data$Température...C. <- NULL
#data2$TempÃ.rature..Â.C. <- NULL
Nous sommes alors à 22 variables.
Mais combien d’observations si l’on retire les lignes avec des manquants ?
which(is.na(data))
data <- na.omit(data)
#which(is.na(data2))
#data2 <- na.omit(data2)
On créé une colonne DateEnJour afin que chaque jour est son numéro
#library(dplyr)
#df %>%
# slice(which.max(as.Date(data$Date)))
data$jour = format(as.Date(data$Date,format="%Y-%m-%d"), format = "%d")
data$mois = ( format(as.Date(data$Date,format="%Y-%m-%d"), format = "%m"))
data$annee = format(as.Date(data$Date,format="%Y-%m-%d"), format = "%Y")
#multiplieur
#c'était un problème de classe de donnée
#data$trente <- replicate(4965, 30)
#data$troisCent <- replicate(4965, 365)
data$jour <- as.numeric(data$jour)
data$mois <- as.numeric(data$mois)
data$annee <- as.numeric(data$annee)
data$DateEnJour <- data$jour + (30 *data$mois) + ( 365*data$annee)
On ne garde que une observation par jour pour un jeu de prévision en effet le nbre d’observation par jour est irrégulier.
prevision <- do.call(rbind, by(data, list(data$DateEnJour),
FUN=function(x) head(x, 1)))
On supprimme les variables créé et la date
data$Date <- NULL
data$DateEnJour <- NULL
data$annee <- NULL
data$mois <- NULL
data$jour <- NULL
prevision$Date <- NULL
prevision$DateEnJour <- NULL
#prevision$annee <- NULL
#prevision$mois <- NULL
prevision$jour <- NULL
On a pas supprimmer le paramètre mois pour faire un test (en utilisant la statistique de student.)
prevision$mois <- ifelse(prevision$mois < 5 ,"Hiver","Pas Hiver")
Il y a t’il une différence significative de moyenne de température entre les relevés en Hiver et les autres ?
#Welch correction inclus dans la fonction t.test()
# test1<-t.test(TempÃ.rature ~ mois, prevision)
# test1
test1<-t.test(Température ~ mois, prevision)
test1##
## Welch Two Sample t-test
##
## data: Température by mois
## t = -26.435, df = 1376.9, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Hiver and group Pas Hiver is not equal to 0
## 95 percent confidence interval:
## -6.995227 -6.028752
## sample estimates:
## mean in group Hiver mean in group Pas Hiver
## 281.355 287.867t est la statistique de student (-26.435) et df est le degré de liberté, la p-value, qui détermine la significativité du test est inférieur à 0.05, on rejette HO : “Les moyennes de températures sont égales” La différence est très significative.
prevision$mois <- NULL
prevision2 <- prevision
prevision$annee <- NULL
Pour le jeu de prévision, nous creons une colonne température mais décaler de 1 jour.
La première méthode utilisé n’était pas correct puisque l’on décalait dans le mauvais sens la colonne avec la méthode lag, puis on supprimait la première observation.
#Mauvaise méthode
#data$TemperatureDecaler <- lag(data$TempÃ.rature)
#on supprimme la première observation
#data = data[-1,]
library(Hmisc)## Le chargement a nécessité le package : lattice## Le chargement a nécessité le package : survival## Le chargement a nécessité le package : Formula## Le chargement a nécessité le package : ggplot2##
## Attachement du package : 'Hmisc'## Les objets suivants sont masqués depuis 'package:base':
##
## format.pval, units#prevision$Precipitation_Jour_Suivant <- Lag(prevision$PrÃ.cipitations.dans.les.24.derniÃ.res.heures, shift = -1)
prevision$Temperature_Jour_Suivant <- Lag(prevision$Température, shift = -1)
#on supprimme la derniere observation
library(dplyr)##
## Attachement du package : 'dplyr'## Les objets suivants sont masqués depuis 'package:Hmisc':
##
## src, summarize## Les objets suivants sont masqués depuis 'package:stats':
##
## filter, lag## Les objets suivants sont masqués depuis 'package:base':
##
## intersect, setdiff, setequal, unionprevision = slice(prevision, 1:(n() -1))
#De même on réduit la taille du nbre d'observation de data
data = slice(data, 1:(n() -662))
dim(prevision)## [1] 1655 22On a donc 1655 observations de 21 variables, cela est suffisant pour la régression multiple.
Bilan des problèmatiques rencontrés avec le jeu de données et leurs résolutions.
Le nombre de variables est très important, les données qualitatives ou manquantes empêche d’effectuer une régression
- On résume les données pour identifier les qualitatives et les supprimmer.
Les données manquantes empêche d’effectuer une régression
- On sélectionne les lignes avec une variables manquantes et on les supprimme.
Si l’on supprimme ces variables, les données seront nulle
- On décide de préalablement supprimmés les variables avec des manquants.
Une fois les variables avec des manquants supprimmés, les données restantes n’expliquent rien, la régression associé a un R² proche de 0
- On regarde combien de manquants sont présents par variables, on classe par ordre décroissant les variables selon leur nombre de manquants, puis on supprimme petit à petit les variables pour avoir le bon équilibre entre variables explicatives et nombre d’observations.
Nous souhaitons réaliser une prévision de la température à la prochaine observation, mais cette variables n’est pas présente
- Création de cette variables par utilisation de la méthode lag.
Le modèle a un R² proche de 0 car les données ne sont pas ordonnées chronologiquement
- On ordonne selon la variable Date
Les prédictions donnent de mauvais résultat, car les observations n’ont pas les même écarts temporelles entre elles, en effet certains jours il y a 0 observation, tandis que d’autres il y a 15.
- On associe à un id à chaque jour en isolant les parties de la variables date puis en effectuant le calcul jour + 30 x mois + 365 x annee, ainsi chaque observations a désormais l’id de son jour. Ensuite on ne garde que la première observation de chaque id, il ne reste plus que une observation par jour.
Statistiques Descriptives
Regardons ce qui caractérise cette base de données, tant au niveau des moyennes, de la variance…
summary(data)## Pression.au.niveau.mer Variation.de.pression.en.3.heures
## Min. : 97170 Min. :-840.000
## 1st Qu.:100850 1st Qu.: -70.000
## Median :101450 Median : 0.000
## Mean :101361 Mean : 7.747
## 3rd Qu.:102000 3rd Qu.: 80.000
## Max. :104370 Max. :1150.000
## Direction.du.vent.moyen.10.mn Vitesse.du.vent.moyen.10.mn Température
## Min. : 0.0 Min. : 0.000 Min. :269.6
## 1st Qu.:180.0 1st Qu.: 2.500 1st Qu.:281.6
## Median :230.0 Median : 3.700 Median :285.1
## Mean :212.8 Mean : 4.014 Mean :285.6
## 3rd Qu.:270.0 3rd Qu.: 5.200 3rd Qu.:289.6
## Max. :360.0 Max. :15.200 Max. :305.4
## Point.de.rosée Humidité Visibilité.horizontale Nebulosité.totale
## Min. :262.9 Min. : 30.00 Min. : 180 Min. : 10.00
## 1st Qu.:278.9 1st Qu.: 72.00 1st Qu.:15000 1st Qu.: 90.00
## Median :282.1 Median : 84.00 Median :30000 Median :100.00
## Mean :282.1 Mean : 80.55 Mean :31057 Mean : 94.53
## 3rd Qu.:285.6 3rd Qu.: 92.00 3rd Qu.:48905 3rd Qu.:100.00
## Max. :294.2 Max. :100.00 Max. :60000 Max. :100.00
## Hauteur.de.la.base.des.nuages.de.l.étage.inférieur Pression.station
## Min. : 30 Min. : 95860
## 1st Qu.: 450 1st Qu.: 99510
## Median : 800 Median :100120
## Mean :1211 Mean :100020
## 3rd Qu.:1250 3rd Qu.:100650
## Max. :7000 Max. :102950
## Variation.de.pression.en.24.heures Température.minimale.du.sol.sur.12.heures
## Min. :-4410.00 Min. :265.1
## 1st Qu.: -435.00 1st Qu.:278.8
## Median : -30.00 Median :282.1
## Mean : -52.56 Mean :282.6
## 3rd Qu.: 350.00 3rd Qu.:286.6
## Max. : 4070.00 Max. :298.9
## Rafale.sur.les.10.dernières.minutes Rafales.sur.une.période
## Min. : 0.000 Min. : 0.600
## 1st Qu.: 4.000 1st Qu.: 4.950
## Median : 6.000 Median : 7.200
## Mean : 6.567 Mean : 7.687
## 3rd Qu.: 8.500 3rd Qu.: 9.800
## Max. :27.400 Max. :27.400
## Précipitations.dans.la.dernière.heure
## Min. :-0.1000
## 1st Qu.: 0.0000
## Median : 0.0000
## Mean : 0.1502
## 3rd Qu.: 0.0000
## Max. :11.4000
## Précipitations.dans.les.3.dernières.heures
## Min. :-0.1000
## 1st Qu.: 0.0000
## Median : 0.0000
## Mean : 0.4124
## 3rd Qu.: 0.2000
## Max. :26.7000
## Précipitations.dans.les.6.dernières.heures
## Min. :-0.1000
## 1st Qu.: 0.0000
## Median : 0.0000
## Mean : 0.8186
## 3rd Qu.: 0.6000
## Max. :38.5000
## Précipitations.dans.les.12.dernières.heures
## Min. :-0.100
## 1st Qu.: 0.000
## Median : 0.200
## Mean : 1.519
## 3rd Qu.: 1.600
## Max. :38.500
## Précipitations.dans.les.24.dernières.heures Hauteur.de.base.2
## Min. :-0.100 Min. : 150
## 1st Qu.: 0.000 1st Qu.: 840
## Median : 0.800 Median :1410
## Mean : 2.862 Mean :2122
## 3rd Qu.: 3.800 3rd Qu.:2580
## Max. :41.200 Max. :8000Nos observations sont complètes, et les variables toutes quantitatives, nous pouvons commencer la régression multiple.
Nous n’avons pas importer en UTF-8 pour économiser de la place, donnons des labels corrects aux colonnes que nous allons utiliser.
#labels <- c(
# "Pression_Mer",
# "Variation_Pression_3h",
# "Direction_vent",
# "Vitesse_vent",
# "Temperature",
# "Point_de_rosee",
# "Humidite",
# "Visibilite_horizontale",
# "Nebulosite_total",
# "Hauteur_Base_Nuage_Inferieur",
# "Pression_Station",
# "Variation_Pression",
# "Temperature_minimale_sol",
# "Rafales",
# "Rafales_total",
# "Precipitation_1h",
# "Precipitation_3h",
# "Precipitation_6h",
# "Precipitation_12h",
# "Precipitation_24h",
# "Hauteur_Base_Nuage_Superieur"
# "Precipitation_Jour_Suivant"
# )
#names(data) <- labels
#label <- c(
# "Pression_Mer",
# "Variation_Pression_3h",
# "Direction_vent",
# "Vitesse_vent",
# "Temperature",
# "Point_de_rosee",
# "Humidite",
# "Visibilite_horizontale",
# "Nebulosite_total",
# "Hauteur_Base_Nuage_Inferieur",
# "Pression_Station",
# "Variation_Pression",
# "Temperature_minimale_sol",
# "Rafales",
# "Rafales_total",
# "Precipitation_1h",
# "Precipitation_3h",
# "Precipitation_6h",
# "Precipitation_12h",
# "Precipitation_24h",
# "Hauteur_Base_Nuage_Superieur",
# "Precipitation_Jour_Suivant"
# )
#names(prevision) <- label
#Avoir toujours les mêmes résultats
# set.seed(1234)
#
# #on sélectionne 20% des lignes
# nb_lignes <- sample(1:nrow(prevision), nrow(prevision)*0.20)
# testing <- prevision[nb_lignes,]
# prevision <- prevision[-nb_lignes,]Notre jeu de donnée pour la regression linéaire est prêt.
## Régression multiple avec toutes les variables explicatives
Nous allons effectuer plusieurs régressions pour essayer d’expliquer la température en fonction des autres variables explicatives.
Nous allons utiliser le R² pour comparer les différentes régressions. Cette mesure indique la proportion de la variance expliquée par le modèle.
- 0 % le modèle n’explique par la variable Y
- 100 % le modèle explique la variabilité de Y lié à la liaison linéaire des variables explicatives entièrement
Tests de regression sur notre jeu complet.
#reg.mul <- lm(Pression.au.niveau.mer~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(Variation.de.pression.en.3.heures~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(Direction.du.vent.moyen.10.mn~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(Vitesse.du.vent.moyen.10.mn~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(Point.de.rosÃ.e~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(HumiditÃ.~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(VisibilitÃ..horizontale~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(NebulositÃ..totale~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(Hauteur.de.la.base.des.nuages.de.l.Ã.tage.infÃ.rieur~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(Pression.station~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(Variation.de.pression.en.24.heures~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(TempÃ.rature.minimale.du.sol.sur.12.heures~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(Rafale.sur.les.10.derniÃ.res.minutes~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(Rafales.sur.une.pÃ.riode~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(PrÃ.cipitations.dans.la.derniÃ.re.heure~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(PrÃ.cipitations.dans.les.3.derniÃ.res.heures~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(PrÃ.cipitations.dans.les.6.derniÃ.res.heures~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(PrÃ.cipitations.dans.les.12.derniÃ.res.heures~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(PrÃ.cipitations.dans.les.24.derniÃ.res.heures~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(Hauteur.de.base.2~.,data)
#shapiro.test(reg.mul$residuals)
#reg.mul <- lm(TempÃ.rature~.,data)
#shapiro.test(reg.mul$residuals)Toutes les régressions testés affichent des résidus non gaussiens. Ceci n’est pas surprenant dans la mesure où nous avons 5000 observations.
Regression avec toutes les variables en testant la normalité des résidus.
reg.mul <- lm(Temperature_Jour_Suivant ~.,prevision)
shapiro.test(reg.mul$residuals)##
## Shapiro-Wilk normality test
##
## data: reg.mul$residuals
## W = 0.99324, p-value = 7.097e-07
On test Ho : “tt les Ɛi suivent N(0,1)” Ici la p-value est inferieur à 5% donc on rejette H0, d’après le test les Ɛi ne suivent pas des lois normales (résidus non gaussiens).
Visualisons les résidus par un histogramme
hist(reg.mul$residuals)
#hist(reg.mul.prevision$residuals)
On réalise l’histogramme des résidus, ils ont la forme d’une Gaussienne.
Donc les résidus semblent suivre une loi gaussienne.
Le test shapiro affiche des résidus non gaussiens, mais ceci peut être expliqué par notre grand nombre d’observations. Cette hypothèse se confirmera lorsqu’on regarde l’histogramme des résidus.
Ainsi nos tests de Fisher et student seront valables.
Quels options sont possibles lorsque notre régression est mauvaise ?
Premièrement, on vérifie par un test de Fisher que le modèle est valide.
Deuxièmement,on peut ajouter la constante au modèle de regression ou d’autres variables explicatives afin de réduire les résidus.
D’autres options sont possibles : une selection de variables (forward ou backward selection), imposer des contraintes au vecteur de régression (cf Lasso ou Ridge), réaliser une régression sur les composantes principales (ACP : variables non correlés)
Continuons avec la régression multiple sur notre jeu de données.
Résultat de la regression
summary(reg.mul)##
## Call:
## lm(formula = Temperature_Jour_Suivant ~ ., data = prevision)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.0470 -2.3822 -0.1275 2.1315 15.1428
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 4.105e+01 4.527e+01
## Pression.au.niveau.mer -2.950e-02 2.976e-02
## Variation.de.pression.en.3.heures -1.155e-03 7.971e-04
## Direction.du.vent.moyen.10.mn -1.708e-03 1.023e-03
## Vitesse.du.vent.moyen.10.mn 1.600e-01 2.104e-01
## Température 2.475e-01 2.186e-01
## Point.de.rosée 1.587e-01 1.925e-01
## Humidité -1.060e-03 4.239e-02
## Visibilité.horizontale 1.004e-05 5.755e-06
## Nebulosité.totale -8.656e-03 8.563e-03
## Hauteur.de.la.base.des.nuages.de.l.étage.inférieur 6.003e-05 1.085e-04
## Pression.station 3.026e-02 3.016e-02
## Variation.de.pression.en.24.heures 4.937e-04 1.665e-04
## Température.minimale.du.sol.sur.12.heures 3.334e-01 3.816e-02
## Rafale.sur.les.10.dernières.minutes -4.767e-02 1.395e-01
## Rafales.sur.une.période -1.467e-01 7.146e-02
## Précipitations.dans.la.dernière.heure 9.053e-02 2.377e-01
## Précipitations.dans.les.3.dernières.heures -1.441e-02 1.957e-01
## Précipitations.dans.les.6.dernières.heures 1.876e-02 1.181e-01
## Précipitations.dans.les.12.dernières.heures -7.490e-03 7.312e-02
## Précipitations.dans.les.24.dernières.heures 8.592e-03 3.331e-02
## Hauteur.de.base.2 2.138e-04 7.763e-05
## t value Pr(>|t|)
## (Intercept) 0.907 0.36468
## Pression.au.niveau.mer -0.991 0.32161
## Variation.de.pression.en.3.heures -1.449 0.14766
## Direction.du.vent.moyen.10.mn -1.670 0.09519 .
## Vitesse.du.vent.moyen.10.mn 0.760 0.44709
## Température 1.132 0.25778
## Point.de.rosée 0.824 0.41000
## Humidité -0.025 0.98006
## Visibilité.horizontale 1.745 0.08114 .
## Nebulosité.totale -1.011 0.31222
## Hauteur.de.la.base.des.nuages.de.l.étage.inférieur 0.553 0.58006
## Pression.station 1.003 0.31587
## Variation.de.pression.en.24.heures 2.966 0.00306 **
## Température.minimale.du.sol.sur.12.heures 8.737 < 2e-16 ***
## Rafale.sur.les.10.dernières.minutes -0.342 0.73266
## Rafales.sur.une.période -2.053 0.04021 *
## Précipitations.dans.la.dernière.heure 0.381 0.70336
## Précipitations.dans.les.3.dernières.heures -0.074 0.94130
## Précipitations.dans.les.6.dernières.heures 0.159 0.87388
## Précipitations.dans.les.12.dernières.heures -0.102 0.91842
## Précipitations.dans.les.24.dernières.heures 0.258 0.79648
## Hauteur.de.base.2 2.754 0.00595 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.515 on 1633 degrees of freedom
## Multiple R-squared: 0.6694, Adjusted R-squared: 0.6651
## F-statistic: 157.5 on 21 and 1633 DF, p-value: < 2.2e-16#summary(reg.mul.prevision)La moyenne des résidus n’est pas affiché car elle est de 0 (ils sont gaussiens).
La colonne Estimate correspond à notre vecteur des Bêta chapeau (notre estimateur de Bêta).
Pour tester leur significativité on regarde les test de student, les t value sont les statistiques de student, si la p-value < 0.05 on rejette H0. Avec H0 : “Bi = 0” > Pour tout i appartenant à (1..21)
Il nous est indiqué que 5 variables sont significatives, il y en a peu car les variables doivent être très corréllés.
Attention, Nous ne pouvons écrire que Temperature du jour suivant = 40 + TempÃ.rature.minimale.du.sol.sur.12.heures * 3.3 + Variation.de.pression.en.24.heures * 5 , il faut refaire le modèle avec uniquement ces variables.
Un R² de 0.67 est trouvé et le R²ajusté = 0.67 également. On test la validité du R² avec un test de fisher (f calculté = 104.6), la p_value est < 0.05 on rejette H0.
Pour rappel : - R² = SCE / SCT - R²ajusté = 1 - ( SCR / (n-p-1) ) / ( SCE / (n-p) )
D’où 57% de la variabilité est expliquée par le modèle de regréssion.
Pour comparer nos régressions,nous utiliserons le R² ajusté.
Réalisons les intervalles de confiance à 95%.
confint(reg.mul)## 2.5 % 97.5 %
## (Intercept) -4.774731e+01 1.298462e+02
## Pression.au.niveau.mer -8.787545e-02 2.886555e-02
## Variation.de.pression.en.3.heures -2.717965e-03 4.087969e-04
## Direction.du.vent.moyen.10.mn -3.715073e-03 2.985685e-04
## Vitesse.du.vent.moyen.10.mn -2.527253e-01 5.727759e-01
## Température -1.812935e-01 6.762108e-01
## Point.de.rosée -2.189583e-01 5.362818e-01
## Humidité -8.420298e-02 8.208397e-02
## Visibilité.horizontale -1.244294e-06 2.132977e-05
## Nebulosité.totale -2.545085e-02 8.138952e-03
## Hauteur.de.la.base.des.nuages.de.l.étage.inférieur -1.527301e-04 2.727909e-04
## Pression.station -2.889892e-02 8.942090e-02
## Variation.de.pression.en.24.heures 1.671774e-04 8.202668e-04
## Température.minimale.du.sol.sur.12.heures 2.585931e-01 4.083004e-01
## Rafale.sur.les.10.dernières.minutes -3.213231e-01 2.259884e-01
## Rafales.sur.une.période -2.868969e-01 -6.558577e-03
## Précipitations.dans.la.dernière.heure -3.757159e-01 5.567828e-01
## Précipitations.dans.les.3.dernières.heures -3.982950e-01 3.694668e-01
## Précipitations.dans.les.6.dernières.heures -2.129717e-01 2.504835e-01
## Précipitations.dans.les.12.dernières.heures -1.509038e-01 1.359237e-01
## Précipitations.dans.les.24.dernières.heures -5.673918e-02 7.392297e-02
## Hauteur.de.base.2 6.153153e-05 3.660643e-04#confint(reg.mul.prevision)Ainsi sur un échantillon de mesures météorologiques de 100 individus, 95 seront dans ces valeurs.
Test d’une Regression sans la Constante :
reg.mulSansC <- lm(Temperature_Jour_Suivant~.-1,prevision)
shapiro.test(reg.mulSansC$residuals)##
## Shapiro-Wilk normality test
##
## data: reg.mulSansC$residuals
## W = 0.9933, p-value = 7.902e-07summary(reg.mulSansC)##
## Call:
## lm(formula = Temperature_Jour_Suivant ~ . - 1, data = prevision)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.0679 -2.3498 -0.1135 2.1267 15.0405
##
## Coefficients:
## Estimate Std. Error
## Pression.au.niveau.mer -3.817e-03 9.111e-03
## Variation.de.pression.en.3.heures -1.161e-03 7.970e-04
## Direction.du.vent.moyen.10.mn -1.712e-03 1.023e-03
## Vitesse.du.vent.moyen.10.mn 1.640e-01 2.104e-01
## Température 3.627e-01 1.778e-01
## Point.de.rosée 1.766e-01 1.915e-01
## Humidité -1.540e-03 4.238e-02
## Visibilité.horizontale 1.024e-05 5.750e-06
## Nebulosité.totale -8.513e-03 8.561e-03
## Hauteur.de.la.base.des.nuages.de.l.étage.inférieur 5.760e-05 1.084e-04
## Pression.station 4.257e-03 9.341e-03
## Variation.de.pression.en.24.heures 4.861e-04 1.663e-04
## Température.minimale.du.sol.sur.12.heures 3.343e-01 3.815e-02
## Rafale.sur.les.10.dernières.minutes -5.016e-02 1.395e-01
## Rafales.sur.une.période -1.432e-01 7.135e-02
## Précipitations.dans.la.dernière.heure 8.498e-02 2.376e-01
## Précipitations.dans.les.3.dernières.heures -5.469e-03 1.955e-01
## Précipitations.dans.les.6.dernières.heures 1.631e-02 1.181e-01
## Précipitations.dans.les.12.dernières.heures -6.704e-03 7.311e-02
## Précipitations.dans.les.24.dernières.heures 1.030e-02 3.325e-02
## Hauteur.de.base.2 2.165e-04 7.757e-05
## t value Pr(>|t|)
## Pression.au.niveau.mer -0.419 0.67532
## Variation.de.pression.en.3.heures -1.457 0.14541
## Direction.du.vent.moyen.10.mn -1.673 0.09443 .
## Vitesse.du.vent.moyen.10.mn 0.779 0.43590
## Température 2.039 0.04157 *
## Point.de.rosée 0.922 0.35659
## Humidité -0.036 0.97103
## Visibilité.horizontale 1.781 0.07512 .
## Nebulosité.totale -0.994 0.32018
## Hauteur.de.la.base.des.nuages.de.l.étage.inférieur 0.531 0.59534
## Pression.station 0.456 0.64864
## Variation.de.pression.en.24.heures 2.924 0.00351 **
## Température.minimale.du.sol.sur.12.heures 8.764 < 2e-16 ***
## Rafale.sur.les.10.dernières.minutes -0.360 0.71918
## Rafales.sur.une.période -2.007 0.04489 *
## Précipitations.dans.la.dernière.heure 0.358 0.72066
## Précipitations.dans.les.3.dernières.heures -0.028 0.97768
## Précipitations.dans.les.6.dernières.heures 0.138 0.89016
## Précipitations.dans.les.12.dernières.heures -0.092 0.92695
## Précipitations.dans.les.24.dernières.heures 0.310 0.75681
## Hauteur.de.base.2 2.791 0.00532 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.515 on 1634 degrees of freedom
## Multiple R-squared: 0.9999, Adjusted R-squared: 0.9998
## F-statistic: 5.212e+05 on 21 and 1634 DF, p-value: < 2.2e-16Le test de shapiro donne des résidus non gaussiens (sans surprise).
On trouve un R² d’environ 1, ce qui est normal.
On refait une régression avec les variables significatives de reg.mul et la constante
reg.mulfin <- lm(Temperature_Jour_Suivant~Variation.de.pression.en.24.heures +
Hauteur.de.base.2 +
Température.minimale.du.sol.sur.12.heures,
# TempÃ.rature.minimale.du.sol.sur.12.heures ,
prevision
)
summary(reg.mulfin)##
## Call:
## lm(formula = Temperature_Jour_Suivant ~ Variation.de.pression.en.24.heures +
## Hauteur.de.base.2 + Température.minimale.du.sol.sur.12.heures,
## data = prevision)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.9386 -2.4950 -0.1546 2.3365 15.6666
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.079e+01 4.887e+00 10.393 < 2e-16
## Variation.de.pression.en.24.heures 7.263e-04 1.456e-04 4.988 6.74e-07
## Hauteur.de.base.2 3.578e-04 4.638e-05 7.714 2.09e-14
## Température.minimale.du.sol.sur.12.heures 8.284e-01 1.733e-02 47.795 < 2e-16
##
## (Intercept) ***
## Variation.de.pression.en.24.heures ***
## Hauteur.de.base.2 ***
## Température.minimale.du.sol.sur.12.heures ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.832 on 1651 degrees of freedom
## Multiple R-squared: 0.6029, Adjusted R-squared: 0.6022
## F-statistic: 835.5 on 3 and 1651 DF, p-value: < 2.2e-16
On trouve un R² de 0.60 et le R² ajusté également.
Touts les paramètres sont significatifs.
shapiro.test(reg.mulfin$residuals)##
## Shapiro-Wilk normality test
##
## data: reg.mulfin$residuals
## W = 0.99068, p-value = 8.87e-09
### Matrice des corrélations des variables explicative
Corrélation entre toutes les variables
#correlation
mcor <- cor(prevision)
#Les Correlations sans visualisation graphique.
mcor
Visualisation graphique des correlations.
#
library(corrplot)## Warning: le package 'corrplot' a été compilé avec la version R 4.1.2## corrplot 0.92 loadedcorrplot(mcor,type="upper",
order="hclust",
tl.col="black",
tl.srt=45,
col = NULL,
tl.pos='n'
)
La matricé était symétrique, on ne garde que le triangle supérieur.
Bleu foncé > corrélation positive forte. Rouge foncé > corrélation négative forte.
On remarque que des groupes entre variables de memes unités se forme tel vitesse du vent avec les données sur les rafales.
Les corrélations les plus importantes sont celles avec Temperature du jour suivant qui est fortement correlées à 3 autres variables.
A noter que nos variables ont des fortes corrélations, nous verrons cet impact avec la régression sur ACP.
On représentation des résidus, on ne doit pas dépasser 5%.
res.m <- rstudent(reg.mulfin)
plot(res.m, pch=20,cex=1, ylab = "Residus", main="",ylim=c(-4,+4))
lines(x=c(-2,0,2),y=c(2,1,2))
Régression pas à pas
Tentons d’obtenir une régression exploitable.
### Ascendant
On commence par une régression à une seule variable
reg1 <- lm(Temperature_Jour_Suivant~1,prevision)
reg1##
## Call:
## lm(formula = Temperature_Jour_Suivant ~ 1, data = prevision)
##
## Coefficients:
## (Intercept)
## 285.8En simplifiant la formule de la régression à une seule variable, il ne reste plus que y(barre), càd la moyenne.
Voyons quelle variable pouvons nous ajouter à cette régression. On effectue un test de Fisher
add1(reg1, scope =~Pression.au.niveau.mer +
Variation.de.pression.en.3.heures +
Direction.du.vent.moyen.10.mn +
Vitesse.du.vent.moyen.10.mn +
# Point.de.rosÃ.e +
# HumiditÃ. +
# VisibilitÃ..horizontale +
# NebulositÃ..totale +
# Hauteur.de.la.base.des.nuages.de.l.Ã.tage.infÃ.rieur+
Point.de.rosée+
Humidité+
Visibilité.horizontale+
Nebulosité.totale+
Hauteur.de.la.base.des.nuages.de.l.étage.inférieur+
Température.minimale.du.sol.sur.12.heures+
Rafale.sur.les.10.dernières.minutes+
Précipitations.dans.les.3.dernières.heures+
Précipitations.dans.les.6.dernières.heures+
Précipitations.dans.les.12.dernières.heures+
Précipitations.dans.les.24.dernières.heures+
Température+
Pression.station +
Variation.de.pression.en.24.heures +
# TempÃ.rature.minimale.du.sol.sur.12.heures +
# Rafale.sur.les.10.derniÃ.res.minutes +
# Rafales.sur.une.pÃ.riode+
# PrÃ.cipitations.dans.la.derniÃ.re.heure +
# PrÃ.cipitations.dans.les.3.derniÃ.res.heures+
# PrÃ.cipitations.dans.les.6.derniÃ.res.heures +
# PrÃ.cipitations.dans.les.12.derniÃ.res.heures +
# PrÃ.cipitations.dans.les.24.derniÃ.res.heures+
Hauteur.de.base.2
# TempÃ.rature
, test="F" )## Single term additions
##
## Model:
## Temperature_Jour_Suivant ~ 1
## Df Sum of Sq RSS AIC
## <none> 61043 5972.9
## Pression.au.niveau.mer 1 338 60705 5965.7
## Variation.de.pression.en.3.heures 1 1 61041 5974.8
## Direction.du.vent.moyen.10.mn 1 294 60749 5966.9
## Vitesse.du.vent.moyen.10.mn 1 990 60053 5947.8
## Point.de.rosée 1 32236 28806 4732.0
## Humidité 1 6429 54613 5790.7
## Visibilité.horizontale 1 3271 57772 5883.7
## Nebulosité.totale 1 1680 59363 5928.7
## Hauteur.de.la.base.des.nuages.de.l.étage.inférieur 1 2032 59011 5918.8
## Température.minimale.du.sol.sur.12.heures 1 35572 25470 4528.3
## Rafale.sur.les.10.dernières.minutes 1 781 60262 5953.6
## Précipitations.dans.les.3.dernières.heures 1 35 61008 5973.9
## Précipitations.dans.les.6.dernières.heures 1 130 60913 5971.3
## Précipitations.dans.les.12.dernières.heures 1 297 60745 5966.8
## Précipitations.dans.les.24.dernières.heures 1 392 60651 5964.2
## Température 1 36690 24353 4454.0
## Pression.station 1 623 60420 5957.9
## Variation.de.pression.en.24.heures 1 33 61009 5974.0
## Hauteur.de.base.2 1 3214 57828 5885.3
## F value Pr(>F)
## <none>
## Pression.au.niveau.mer 9.1989 0.002459 **
## Variation.de.pression.en.3.heures 0.0342 0.853387
## Direction.du.vent.moyen.10.mn 7.9983 0.004738 **
## Vitesse.du.vent.moyen.10.mn 27.2496 2.014e-07 ***
## Point.de.rosée 1849.8304 < 2.2e-16 ***
## Humidité 194.5943 < 2.2e-16 ***
## Visibilité.horizontale 93.5965 < 2.2e-16 ***
## Nebulosité.totale 46.7708 1.119e-11 ***
## Hauteur.de.la.base.des.nuages.de.l.étage.inférieur 56.9183 7.459e-14 ***
## Température.minimale.du.sol.sur.12.heures 2308.6311 < 2.2e-16 ***
## Rafale.sur.les.10.dernières.minutes 21.4153 3.986e-06 ***
## Précipitations.dans.les.3.dernières.heures 0.9419 0.331942
## Précipitations.dans.les.6.dernières.heures 3.5283 0.060507 .
## Précipitations.dans.les.12.dernières.heures 8.0947 0.004494 **
## Précipitations.dans.les.24.dernières.heures 10.6823 0.001104 **
## Température 2490.4048 < 2.2e-16 ***
## Pression.station 17.0354 3.852e-05 ***
## Variation.de.pression.en.24.heures 0.8996 0.343017
## Hauteur.de.base.2 91.8777 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1On remarque avec les différentes p-values que de nombreuses variables sont éligibles à la rentrée dans le modèle.
La F value nous indique si l’augmentation du R² est significative. Par exemple, pour la variable “Pression.au.niveau.mer”, l’hypothèse H0 est “Bêta^Pression.au.niveau.mer = 0” On rejette cette hypothèse si la p-value est inférieur à 0.05, Donc Bêta^Pression.au.niveau.mer est significatif.
Si cela est significatif on peut introduire la variable dans le modèle.
On remarque que plusieurs variables affiche une pvalue de 2.2e-16 , càd le nbre le plus bas affichable par R.
On départage ces variables en prenant en compte la F value.
On prend la f-value la plus grande,, ce n’est pas préciser mais il semble qu’il s’agisse du coefficient corrigé (lorsque n/k > 40) cela permet de maximiser la qualité de l’ajustement en pénalisant l’over-fitting.
résumé over fitting : très peu voir aucune erreurs sur notre jeu d’entrainement mais beaucoup d’erreurs avec un jeu de test.
Température a la F value la plus forte, sa pvalue permet de rejeter H0 On créé un nouveau modèle avec cette variable ajouter.
On introduit les variables jusqu’à qu’aucune n’ai une augmentation du R² significative lors de son introduction dans le modèle.
add1(update(reg1,~.+Température
+Température.minimale.du.sol.sur.12.heures
+Rafale.sur.les.10.dernières.minutes
+Hauteur.de.base.2
+Pression.station
+Point.de.rosée
+Variation.de.pression.en.24.heures
+Rafales.sur.une.période),
scope =~Pression.au.niveau.mer +
Variation.de.pression.en.3.heures +
Direction.du.vent.moyen.10.mn +
Vitesse.du.vent.moyen.10.mn +
# Point.de.rosÃ.e +
# HumiditÃ. +
# VisibilitÃ..horizontale +
# NebulositÃ..totale +
# Hauteur.de.la.base.des.nuages.de.l.Ã.tage.infÃ.rieur+
Point.de.rosée+
Humidité+
Visibilité.horizontale+
Nebulosité.totale+
Hauteur.de.la.base.des.nuages.de.l.étage.inférieur+
Température.minimale.du.sol.sur.12.heures+
Rafale.sur.les.10.dernières.minutes+
Précipitations.dans.les.3.dernières.heures+
Précipitations.dans.les.6.dernières.heures+
Précipitations.dans.les.12.dernières.heures+
Précipitations.dans.les.24.dernières.heures+
Température+
Rafales.sur.une.période+
Pression.station +
Variation.de.pression.en.24.heures +
# TempÃ.rature.minimale.du.sol.sur.12.heures +
# Rafale.sur.les.10.derniÃ.res.minutes +
# Rafales.sur.une.pÃ.riode+
# PrÃ.cipitations.dans.la.derniÃ.re.heure +
# PrÃ.cipitations.dans.les.3.derniÃ.res.heures+
# PrÃ.cipitations.dans.les.6.derniÃ.res.heures +
# PrÃ.cipitations.dans.les.12.derniÃ.res.heures +
# PrÃ.cipitations.dans.les.24.derniÃ.res.heures+
Hauteur.de.base.2
# TempÃ.rature
, test="F" )add1(update(reg1,~.+Température
+Température.minimale.du.sol.sur.12.heures
+Rafale.sur.les.10.dernières.minutes
+Hauteur.de.base.2
+Pression.station
+Point.de.rosée
+Variation.de.pression.en.24.heures
+Rafales.sur.une.période
+Variation.de.pression.en.3.heures),
scope =~Pression.au.niveau.mer +
Variation.de.pression.en.3.heures +
Direction.du.vent.moyen.10.mn +
Vitesse.du.vent.moyen.10.mn +
# Point.de.rosÃ.e +
# HumiditÃ. +
# VisibilitÃ..horizontale +
# NebulositÃ..totale +
# Hauteur.de.la.base.des.nuages.de.l.Ã.tage.infÃ.rieur+
Point.de.rosée+
Humidité+
Visibilité.horizontale+
Nebulosité.totale+
Hauteur.de.la.base.des.nuages.de.l.étage.inférieur+
Température.minimale.du.sol.sur.12.heures+
Rafale.sur.les.10.dernières.minutes+
Précipitations.dans.les.3.dernières.heures+
Précipitations.dans.les.6.dernières.heures+
Précipitations.dans.les.12.dernières.heures+
Précipitations.dans.les.24.dernières.heures+
Température+
Rafales.sur.une.période+
Pression.station +
Variation.de.pression.en.24.heures +
# TempÃ.rature.minimale.du.sol.sur.12.heures +
# Rafale.sur.les.10.derniÃ.res.minutes +
# Rafales.sur.une.pÃ.riode+
# PrÃ.cipitations.dans.la.derniÃ.re.heure +
# PrÃ.cipitations.dans.les.3.derniÃ.res.heures+
# PrÃ.cipitations.dans.les.6.derniÃ.res.heures +
# PrÃ.cipitations.dans.les.12.derniÃ.res.heures +
# PrÃ.cipitations.dans.les.24.derniÃ.res.heures+
Hauteur.de.base.2
# TempÃ.rature
, test="F" )Aucune variable à introduire, on arrête notre forward selection ici.
Modele de la forward regression
reg.mul.forward <- lm(Temperature_Jour_Suivant~Température
+Température.minimale.du.sol.sur.12.heures
+Rafale.sur.les.10.dernières.minutes
+Hauteur.de.base.2
+Pression.station
+Point.de.rosée
+Variation.de.pression.en.24.heures
+Rafales.sur.une.période
+Variation.de.pression.en.3.heures
,
prevision
)
summary(reg.mul.forward)##
## Call:
## lm(formula = Temperature_Jour_Suivant ~ Température + Température.minimale.du.sol.sur.12.heures +
## Rafale.sur.les.10.dernières.minutes + Hauteur.de.base.2 +
## Pression.station + Point.de.rosée + Variation.de.pression.en.24.heures +
## Rafales.sur.une.période + Variation.de.pression.en.3.heures,
## data = prevision)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.1638 -2.3741 -0.1582 2.1165 15.3027
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 2.454e+00 1.265e+01 0.194
## Température 4.070e-01 3.100e-02 13.129
## Température.minimale.du.sol.sur.12.heures 3.348e-01 3.771e-02 8.877
## Rafale.sur.les.10.dernières.minutes 2.627e-02 7.326e-02 0.359
## Hauteur.de.base.2 2.798e-04 4.671e-05 5.990
## Pression.station 3.040e-04 1.145e-04 2.655
## Point.de.rosée 1.491e-01 4.281e-02 3.483
## Variation.de.pression.en.24.heures 5.048e-04 1.600e-04 3.154
## Rafales.sur.une.période -1.259e-01 6.882e-02 -1.829
## Variation.de.pression.en.3.heures -1.350e-03 7.617e-04 -1.772
## Pr(>|t|)
## (Intercept) 0.846270
## Température < 2e-16 ***
## Température.minimale.du.sol.sur.12.heures < 2e-16 ***
## Rafale.sur.les.10.dernières.minutes 0.719921
## Hauteur.de.base.2 2.57e-09 ***
## Pression.station 0.008015 **
## Point.de.rosée 0.000509 ***
## Variation.de.pression.en.24.heures 0.001639 **
## Rafales.sur.une.période 0.067531 .
## Variation.de.pression.en.3.heures 0.076598 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.511 on 1645 degrees of freedom
## Multiple R-squared: 0.6677, Adjusted R-squared: 0.6659
## F-statistic: 367.3 on 9 and 1645 DF, p-value: < 2.2e-16
Avec cette méthode on peut écrire l’équation du modèle. Y : la température et Bêta chapeau le vecteur issues de la colonne Estimate de dimension 9 (avec la constante).
On a un R² de 0.67 et un R² également.
On test les résidus
shapiro.test(reg.mul.forward$residuals)##
## Shapiro-Wilk normality test
##
## data: reg.mul.forward$residuals
## W = 0.99335, p-value = 8.641e-07
Les résidus ne sont toujours pas Gaussiens par le test de Shapiro.
Méthode forwards (procédure automatique).
full <- lm(Temperature_Jour_Suivant ~.,prevision)
null <- lm(Temperature_Jour_Suivant ~1,prevision)
forw <- step(null,scope=list(lower=null, upper = full),direction = "forward",trace=0)
formula(forw)## Temperature_Jour_Suivant ~ Température + Température.minimale.du.sol.sur.12.heures +
## Rafales.sur.une.période + Hauteur.de.base.2 + Pression.station +
## Point.de.rosée + Variation.de.pression.en.24.heures + Variation.de.pression.en.3.heures +
## Visibilité.horizontale + Direction.du.vent.moyen.10.mnsummary(lm(forw,prevision))##
## Call:
## lm(formula = forw, data = prevision)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.1651 -2.3510 -0.1287 2.1048 15.1107
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -8.667e-01 1.273e+01 -0.068
## Température 3.875e-01 3.355e-02 11.551
## Température.minimale.du.sol.sur.12.heures 3.341e-01 3.775e-02 8.850
## Rafales.sur.une.période -1.014e-01 2.786e-02 -3.640
## Hauteur.de.base.2 2.647e-04 4.753e-05 5.569
## Pression.station 3.350e-04 1.154e-04 2.903
## Point.de.rosée 1.705e-01 4.478e-02 3.808
## Variation.de.pression.en.24.heures 5.012e-04 1.600e-04 3.133
## Variation.de.pression.en.3.heures -1.214e-03 7.682e-04 -1.581
## Visibilité.horizontale 9.216e-06 5.364e-06 1.718
## Direction.du.vent.moyen.10.mn -1.582e-03 1.010e-03 -1.566
## Pr(>|t|)
## (Intercept) 0.945710
## Température < 2e-16 ***
## Température.minimale.du.sol.sur.12.heures < 2e-16 ***
## Rafales.sur.une.période 0.000281 ***
## Hauteur.de.base.2 2.98e-08 ***
## Pression.station 0.003742 **
## Point.de.rosée 0.000145 ***
## Variation.de.pression.en.24.heures 0.001758 **
## Variation.de.pression.en.3.heures 0.114178
## Visibilité.horizontale 0.085936 .
## Direction.du.vent.moyen.10.mn 0.117463
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.507 on 1644 degrees of freedom
## Multiple R-squared: 0.6687, Adjusted R-squared: 0.6667
## F-statistic: 331.8 on 10 and 1644 DF, p-value: < 2.2e-16shapiro.test(lm(forw,prevision)$residuals)##
## Shapiro-Wilk normality test
##
## data: lm(forw, prevision)$residuals
## W = 0.99335, p-value = 8.646e-07
Entre la procédure automatique et manuelle on a quelques variables différentes, mais on trouve un R² similaire.
### Descendant Méthode Backward selection (procédure manuelle).
On commence par voir quel paramètre on retirera en premier.
On utilise le test de Fisher pour cela pour déterminer celui qui est la moins significative > celle qui fait diminuer le R².
On rappel que pour le test de fisher, la statistique calculée est : - ( SCE / n ) / ( SCR / (n-p-1))
On regarde la variable qui a la plus petite F Value, et on la supprime du modèle.
lm.full <- lm(Temperature_Jour_Suivant ~.,prevision)
drop1(lm.full,test="F")## Single term deletions
##
## Model:
## Temperature_Jour_Suivant ~ Pression.au.niveau.mer + Variation.de.pression.en.3.heures +
## Direction.du.vent.moyen.10.mn + Vitesse.du.vent.moyen.10.mn +
## Température + Point.de.rosée + Humidité + Visibilité.horizontale +
## Nebulosité.totale + Hauteur.de.la.base.des.nuages.de.l.étage.inférieur +
## Pression.station + Variation.de.pression.en.24.heures + Température.minimale.du.sol.sur.12.heures +
## Rafale.sur.les.10.dernières.minutes + Rafales.sur.une.période +
## Précipitations.dans.la.dernière.heure + Précipitations.dans.les.3.dernières.heures +
## Précipitations.dans.les.6.dernières.heures + Précipitations.dans.les.12.dernières.heures +
## Précipitations.dans.les.24.dernières.heures + Hauteur.de.base.2
## Df Sum of Sq RSS AIC
## <none> 20181 4183.0
## Pression.au.niveau.mer 1 12.15 20193 4182.0
## Variation.de.pression.en.3.heures 1 25.93 20207 4183.2
## Direction.du.vent.moyen.10.mn 1 34.45 20215 4183.9
## Vitesse.du.vent.moyen.10.mn 1 7.15 20188 4181.6
## Température 1 15.84 20197 4182.3
## Point.de.rosée 1 8.39 20189 4181.7
## Humidité 1 0.01 20181 4181.0
## Visibilité.horizontale 1 37.64 20218 4184.1
## Nebulosité.totale 1 12.63 20193 4182.1
## Hauteur.de.la.base.des.nuages.de.l.étage.inférieur 1 3.78 20185 4181.4
## Pression.station 1 12.44 20193 4182.1
## Variation.de.pression.en.24.heures 1 108.69 20290 4189.9
## Température.minimale.du.sol.sur.12.heures 1 943.45 21124 4256.7
## Rafale.sur.les.10.dernières.minutes 1 1.44 20182 4181.2
## Rafales.sur.une.période 1 52.10 20233 4185.3
## Précipitations.dans.la.dernière.heure 1 1.79 20183 4181.2
## Précipitations.dans.les.3.dernières.heures 1 0.07 20181 4181.0
## Précipitations.dans.les.6.dernières.heures 1 0.31 20181 4181.1
## Précipitations.dans.les.12.dernières.heures 1 0.13 20181 4181.1
## Précipitations.dans.les.24.dernières.heures 1 0.82 20182 4181.1
## Hauteur.de.base.2 1 93.73 20275 4188.7
## F value Pr(>F)
## <none>
## Pression.au.niveau.mer 0.9830 0.321612
## Variation.de.pression.en.3.heures 2.0983 0.147658
## Direction.du.vent.moyen.10.mn 2.7876 0.095189 .
## Vitesse.du.vent.moyen.10.mn 0.5783 0.447095
## Température 1.2815 0.257779
## Point.de.rosée 0.6792 0.409995
## Humidité 0.0006 0.980062
## Visibilité.horizontale 3.0457 0.081140 .
## Nebulosité.totale 1.0219 0.312215
## Hauteur.de.la.base.des.nuages.de.l.étage.inférieur 0.3063 0.580056
## Pression.station 1.0066 0.315870
## Variation.de.pression.en.24.heures 8.7947 0.003065 **
## Température.minimale.du.sol.sur.12.heures 76.3426 < 2.2e-16 ***
## Rafale.sur.les.10.dernières.minutes 0.1167 0.732655
## Rafales.sur.une.période 4.2156 0.040213 *
## Précipitations.dans.la.dernière.heure 0.1451 0.703360
## Précipitations.dans.les.3.dernières.heures 0.0054 0.941299
## Précipitations.dans.les.6.dernières.heures 0.0252 0.873881
## Précipitations.dans.les.12.dernières.heures 0.0105 0.918421
## Précipitations.dans.les.24.dernières.heures 0.0665 0.796476
## Hauteur.de.base.2 7.5847 0.005952 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1à chaque itération, on retire la p-value de F la plus faible, ici Humidité
Toutes les variables sont significatives, on a le modele backward selection complet.
summary(update(lm.full,~.-Vitesse.du.vent.moyen.10.mn
-Température-Hauteur.de.base.2
-Variation.de.pression.en.3.heures
-Précipitations.dans.les.24.dernières.heures
-Direction.du.vent.moyen.10.mn
-Précipitations.dans.les.3.dernières.heures
-Rafale.sur.les.10.dernières.minutes
-Visibilité.horizontale
-Pression.station
-Précipitations.dans.les.6.dernières.heures
-Précipitations.dans.les.12.dernières.heures
-Pression.au.niveau.mer
-Précipitations.dans.la.dernière.heure ))##
## Call:
## lm(formula = Temperature_Jour_Suivant ~ Point.de.rosée + Humidité +
## Nebulosité.totale + Hauteur.de.la.base.des.nuages.de.l.étage.inférieur +
## Variation.de.pression.en.24.heures + Température.minimale.du.sol.sur.12.heures +
## Rafales.sur.une.période, data = prevision)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.9053 -2.3955 -0.1457 2.1538 15.1599
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 4.121e+01 5.290e+00
## Point.de.rosée 5.823e-01 4.141e-02
## Humidité -9.314e-02 7.071e-03
## Nebulosité.totale -1.913e-02 7.485e-03
## Hauteur.de.la.base.des.nuages.de.l.étage.inférieur 3.037e-04 7.221e-05
## Variation.de.pression.en.24.heures 4.872e-04 1.393e-04
## Température.minimale.du.sol.sur.12.heures 3.196e-01 3.740e-02
## Rafales.sur.une.période -1.562e-01 2.664e-02
## t value Pr(>|t|)
## (Intercept) 7.790 1.18e-14 ***
## Point.de.rosée 14.061 < 2e-16 ***
## Humidité -13.172 < 2e-16 ***
## Nebulosité.totale -2.556 0.010674 *
## Hauteur.de.la.base.des.nuages.de.l.étage.inférieur 4.206 2.74e-05 ***
## Variation.de.pression.en.24.heures 3.498 0.000482 ***
## Température.minimale.du.sol.sur.12.heures 8.545 < 2e-16 ***
## Rafales.sur.une.période -5.866 5.39e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.532 on 1647 degrees of freedom
## Multiple R-squared: 0.6633, Adjusted R-squared: 0.6619
## F-statistic: 463.6 on 7 and 1647 DF, p-value: < 2.2e-16shapiro.test(update(lm.full,~.-Vitesse.du.vent.moyen.10.mn
-Température-Hauteur.de.base.2
-Variation.de.pression.en.3.heures
-Précipitations.dans.les.24.dernières.heures
-Direction.du.vent.moyen.10.mn
-Précipitations.dans.les.3.dernières.heures
-Rafale.sur.les.10.dernières.minutes
-Visibilité.horizontale
-Pression.station
-Précipitations.dans.les.6.dernières.heures
-Précipitations.dans.les.12.dernières.heures
-Pression.au.niveau.mer
-Précipitations.dans.la.dernière.heure )$residuals)##
## Shapiro-Wilk normality test
##
## data: update(lm.full, ~. - Vitesse.du.vent.moyen.10.mn - Température - Hauteur.de.base.2 - Variation.de.pression.en.3.heures - Précipitations.dans.les.24.dernières.heures - Direction.du.vent.moyen.10.mn - Précipitations.dans.les.3.dernières.heures - Rafale.sur.les.10.dernières.minutes - Visibilité.horizontale - Pression.station - Précipitations.dans.les.6.dernières.heures - Précipitations.dans.les.12.dernières.heures - Pression.au.niveau.mer - Précipitations.dans.la.dernière.heure)$residuals
## W = 0.99179, p-value = 5.421e-08# summary(update(lm.full,~.-HumiditÃ.
# -PrÃ.cipitations.dans.les.3.derniÃ.res.heures
# -PrÃ.cipitations.dans.les.12.derniÃ.res.heures
# -PrÃ.cipitations.dans.les.6.derniÃ.res.heures
# -Rafale.sur.les.10.derniÃ.res.minutes
# -PrÃ.cipitations.dans.les.24.derniÃ.res.heures
# -Hauteur.de.la.base.des.nuages.de.l.Ã.tage.infÃ.rieur
# -PrÃ.cipitations.dans.la.derniÃ.re.heure
# -Vitesse.du.vent.moyen.10.mn
# -NebulositÃ..totale
# -Pression.au.niveau.mer
# -Direction.du.vent.moyen.10.mn
# -VisibilitÃ..horizontale ))
#
#
# shapiro.test(update(lm.full,~.-HumiditÃ.
# -PrÃ.cipitations.dans.les.3.derniÃ.res.heures
# -PrÃ.cipitations.dans.les.12.derniÃ.res.heures
# -PrÃ.cipitations.dans.les.6.derniÃ.res.heures
# -Rafale.sur.les.10.derniÃ.res.minutes
# -PrÃ.cipitations.dans.les.24.derniÃ.res.heures
# -Hauteur.de.la.base.des.nuages.de.l.Ã.tage.infÃ.rieur
# -PrÃ.cipitations.dans.la.derniÃ.re.heure
# -Vitesse.du.vent.moyen.10.mn
# -NebulositÃ..totale
# -Pression.au.niveau.mer
# -Direction.du.vent.moyen.10.mn
# -VisibilitÃ..horizontale )$residuals)Avec cette méthode on a un R² de 0.67, c’est également le cas du R² ajusté.
On peut écrire le modèle avec le Béta chapeau : la colonne estimate.
Test par réalisation automatique de la backward selection.
back <- step(full, direction="backward")formula(back)## Temperature_Jour_Suivant ~ Variation.de.pression.en.3.heures +
## Direction.du.vent.moyen.10.mn + Température + Point.de.rosée +
## Visibilité.horizontale + Pression.station + Variation.de.pression.en.24.heures +
## Température.minimale.du.sol.sur.12.heures + Rafales.sur.une.période +
## Hauteur.de.base.2reg.mul.backward <- lm(back,prevision)
summary(reg.mul.backward)##
## Call:
## lm(formula = back, data = prevision)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.1651 -2.3510 -0.1287 2.1048 15.1107
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -8.667e-01 1.273e+01 -0.068
## Variation.de.pression.en.3.heures -1.214e-03 7.682e-04 -1.581
## Direction.du.vent.moyen.10.mn -1.582e-03 1.010e-03 -1.566
## Température 3.875e-01 3.355e-02 11.551
## Point.de.rosée 1.705e-01 4.478e-02 3.808
## Visibilité.horizontale 9.216e-06 5.364e-06 1.718
## Pression.station 3.350e-04 1.154e-04 2.903
## Variation.de.pression.en.24.heures 5.012e-04 1.600e-04 3.133
## Température.minimale.du.sol.sur.12.heures 3.341e-01 3.775e-02 8.850
## Rafales.sur.une.période -1.014e-01 2.786e-02 -3.640
## Hauteur.de.base.2 2.647e-04 4.753e-05 5.569
## Pr(>|t|)
## (Intercept) 0.945710
## Variation.de.pression.en.3.heures 0.114178
## Direction.du.vent.moyen.10.mn 0.117463
## Température < 2e-16 ***
## Point.de.rosée 0.000145 ***
## Visibilité.horizontale 0.085936 .
## Pression.station 0.003742 **
## Variation.de.pression.en.24.heures 0.001758 **
## Température.minimale.du.sol.sur.12.heures < 2e-16 ***
## Rafales.sur.une.période 0.000281 ***
## Hauteur.de.base.2 2.98e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.507 on 1644 degrees of freedom
## Multiple R-squared: 0.6687, Adjusted R-squared: 0.6667
## F-statistic: 331.8 on 10 and 1644 DF, p-value: < 2.2e-16shapiro.test(reg.mul.backward$residuals)##
## Shapiro-Wilk normality test
##
## data: reg.mul.backward$residuals
## W = 0.99335, p-value = 8.646e-07Les variables sélectionnés sont en grande partie semblables et donne un R² similaire.
Régression sur ACP
Continuons à rechercher une régression valide.
On réalise une ACP normée avec 21 composantes principales.
Les valeurs sont centrées et réduites : scale.unit = True, Nombre de composantes principales 21 : ncp.
library(FactoMineR)
res.pca = PCA(prevision[,1:21],scale.unit=TRUE, ncp=21,graph=T)## Warning: ggrepel: 5 unlabeled data points (too many overlaps). Consider
## increasing max.overlapsplot.PCA(res.pca, axes=c(1, 2), choix="ind", habillage=5,label="var",graph.type = "ggplot")
On sort les valeurs propre de la matrice de corrélation
res.pca$eig## eigenvalue percentage of variance cumulative percentage of variance
## comp 1 4.975185e+00 2.369136e+01 23.69136
## comp 2 3.185857e+00 1.517075e+01 38.86211
## comp 3 2.833054e+00 1.349073e+01 52.35284
## comp 4 2.079102e+00 9.900487e+00 62.25333
## comp 5 1.526669e+00 7.269851e+00 69.52318
## comp 6 1.257222e+00 5.986772e+00 75.50995
## comp 7 1.000499e+00 4.764281e+00 80.27423
## comp 8 8.150179e-01 3.881038e+00 84.15527
## comp 9 7.399857e-01 3.523741e+00 87.67901
## comp 10 7.286337e-01 3.469684e+00 91.14869
## comp 11 4.989004e-01 2.375716e+00 93.52441
## comp 12 4.585494e-01 2.183568e+00 95.70798
## comp 13 2.877600e-01 1.370286e+00 97.07826
## comp 14 1.724224e-01 8.210592e-01 97.89932
## comp 15 1.447523e-01 6.892966e-01 98.58862
## comp 16 1.139587e-01 5.426604e-01 99.13128
## comp 17 8.263508e-02 3.935004e-01 99.52478
## comp 18 7.639371e-02 3.637796e-01 99.88856
## comp 19 2.032383e-02 9.678014e-02 99.98534
## comp 20 3.073026e-03 1.463346e-02 99.99997
## comp 21 5.351784e-06 2.548468e-05 100.00000On remarque que avec les pts projeté sur la composante 1, on récupère 23% de la variance total (lamba / ncp) avec eigenvalue la valeur de lambda.
On détermine la matrice des scores, toutes les dims sont orthogonaux.
res.pca$ind$coord
novariables<-res.pca$ind$coord
previsionACP <- as.data.frame(res.pca$ind$coord)
previsionACP$Temperature_Jour_Suivant <- prevision$Temperature_Jour_Suivant
Méthode de régression sur ACP :
Méthode Forward :
Sur toutes les données, et sur la constante.
lm.full <- lm(Temperature_Jour_Suivant ~.,previsionACP)
lm.full##
## Call:
## lm(formula = Temperature_Jour_Suivant ~ ., data = previsionACP)
##
## Coefficients:
## (Intercept) Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## 285.82577 -0.50683 2.18175 1.55375 0.50840 -0.32164
## Dim.6 Dim.7 Dim.8 Dim.9 Dim.10 Dim.11
## -0.36316 -0.47258 0.41756 0.12695 -0.18470 0.42103
## Dim.12 Dim.13 Dim.14 Dim.15 Dim.16 Dim.17
## 0.09681 -0.06155 0.19354 -0.03611 -0.27072 0.46745
## Dim.18 Dim.19 Dim.20 Dim.21
## -0.25438 -0.31357 1.01701 -37.20705lm.null <- lm(Temperature_Jour_Suivant ~1,previsionACP)
lm.null##
## Call:
## lm(formula = Temperature_Jour_Suivant ~ 1, data = previsionACP)
##
## Coefficients:
## (Intercept)
## 285.8
Le principe est de tester pour chaque dim trouvé, si elle est significative par un test de fisher.
On ajoute à notre modèle celui qui est la plus significative.
Dans la première itération : * -> Dim 2 à une F value de 827.3598 pour une probabilité d’être supérieur à F < 2.2e-16. (donc très significative) * -> On l’ajoute au modèle * -> on passe à l’itération suivante, jusqu’à ne plus avoir de variables significatives.
En refaisant cette regression en prenant soit le plus petit RSS ou soit les dimensions dans l’ordre, le résultat est le meme.
On a donc :
reg.ACP <- lm(Temperature_Jour_Suivant~Dim.1 + Dim.2 + Dim.3 + Dim.4 + Dim.5 + Dim.6 + Dim.7 + Dim.8 + Dim.11 + Dim.10,previsionACP)
summary(reg.ACP)##
## Call:
## lm(formula = Temperature_Jour_Suivant ~ Dim.1 + Dim.2 + Dim.3 +
## Dim.4 + Dim.5 + Dim.6 + Dim.7 + Dim.8 + Dim.11 + Dim.10,
## data = previsionACP)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.0378 -2.3754 -0.1332 2.0948 15.2685
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 285.82577 0.08636 3309.607 < 2e-16 ***
## Dim.1 -0.50683 0.03872 -13.090 < 2e-16 ***
## Dim.2 2.18175 0.04839 45.091 < 2e-16 ***
## Dim.3 1.55375 0.05131 30.282 < 2e-16 ***
## Dim.4 0.50840 0.05989 8.488 < 2e-16 ***
## Dim.5 -0.32164 0.06990 -4.602 4.51e-06 ***
## Dim.6 -0.36316 0.07702 -4.715 2.62e-06 ***
## Dim.7 -0.47258 0.08634 -5.473 5.10e-08 ***
## Dim.8 0.41756 0.09566 4.365 1.35e-05 ***
## Dim.11 0.42103 0.12227 3.443 0.000589 ***
## Dim.10 -0.18470 0.10117 -1.826 0.068090 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.513 on 1644 degrees of freedom
## Multiple R-squared: 0.6676, Adjusted R-squared: 0.6655
## F-statistic: 330.1 on 10 and 1644 DF, p-value: < 2.2e-16
On a donc un modèle avec l’équation Y (température) et le Bêta chapeau qui est le vecteur estimate.
On obtient un R² de 0.67 et un R² de même valeur.
shapiro.test(reg.ACP$residuals)##
## Shapiro-Wilk normality test
##
## data: reg.ACP$residuals
## W = 0.99281, p-value = 3.202e-07
## Régression Ridge
Quel est la différence avec la régression que nous venons d’effectuer ?
la MLR estime les paramètres sans biais. Tandis que Lasso applique une contrainte linéaire et Ridge une contrainte non linéaire à ces paramètres à estimer.
Cependant, en MLR les coefficients peuvent être très sensibles à l’échantillonnage, et des coefficients censés être positifs pour avoir un sens (ex: une surface) sont négatifs.
Avec les méthodes Lasso et Ridge, en posant des contraintes, on évite ces problèmes.
library(MASS)##
## Attachement du package : 'MASS'## L'objet suivant est masqué depuis 'package:dplyr':
##
## selectmod.ridge = lm.ridge(Temperature_Jour_Suivant ~.,prevision, lambda=seq(0,200,0.01))
plot(mod.ridge)
legend("topright",legend=rownames(mod.ridge$coef),col=1:21,lty=1:3)
On optimise la régression Ridge
select(mod.ridge)## modified HKB estimator is 0.1683841
## modified L-W estimator is 9.510109
## smallest value of GCV at 101.01
Création du modèle optimiser
mod.ridgeopt =lm.ridge(Temperature_Jour_Suivant ~. , data=prevision, lambda = 69.36 )
Coefficient du modele
coeff = coef(mod.ridgeopt)
coeff##
## 4.250563e+00
## Pression.au.niveau.mer
## 1.390285e-04
## Variation.de.pression.en.3.heures
## -1.124248e-03
## Direction.du.vent.moyen.10.mn
## -1.484909e-03
## Vitesse.du.vent.moyen.10.mn
## 4.987332e-02
## Température
## 2.570599e-01
## Point.de.rosée
## 3.075717e-01
## Humidité
## -2.864446e-02
## Visibilité.horizontale
## 9.570773e-06
## Nebulosité.totale
## -9.622582e-03
## Hauteur.de.la.base.des.nuages.de.l.étage.inférieur
## 7.692027e-05
## Pression.station
## 2.092797e-04
## Variation.de.pression.en.24.heures
## 4.536590e-04
## Température.minimale.du.sol.sur.12.heures
## 3.171223e-01
## Rafale.sur.les.10.dernières.minutes
## -2.173624e-02
## Rafales.sur.une.période
## -1.129750e-01
## Précipitations.dans.la.dernière.heure
## 5.559390e-02
## Précipitations.dans.les.3.dernières.heures
## 1.288300e-02
## Précipitations.dans.les.6.dernières.heures
## 1.366126e-02
## Précipitations.dans.les.12.dernières.heures
## -6.823482e-03
## Précipitations.dans.les.24.dernières.heures
## 6.884081e-03
## Hauteur.de.base.2
## 2.047509e-04
Test du modèle obtenu
summary(mod.ridgeopt)## Length Class Mode
## coef 21 -none- numeric
## scales 21 -none- numeric
## Inter 1 -none- numeric
## lambda 1 -none- numeric
## ym 1 -none- numeric
## xm 21 -none- numeric
## GCV 1 -none- numeric
## kHKB 1 -none- numeric
## kLW 1 -none- numeric#shapiro.test(mod.ridgeopt$residuals)
Utilisation d’une autre librairie pour réaliser une prédiction de la régression Ridge.
Cette librairie ne fonctionne pas avec des data frame, on convertit nos données
var <- prevision[,names(prevision) != 'Temperature_Jour_Suivant']
var <- as.matrix(var)
On utilise la validation croisé pour éviter le sur-apprentissage.
library(glmnet)## Warning: le package 'glmnet' a été compilé avec la version R 4.1.2## Le chargement a nécessité le package : Matrix## Loaded glmnet 4.1-3ridge_opti = cv.glmnet(x = var,
y = prevision$Temperature_Jour_Suivant,
family = 'gaussian',
nfolds = 10
)
ridge_lambda <- ridge_opti$lambda.min
ridge_lambda
summary(ridge_opti)
On réalise une prédiction sur nos variables à partir du modèle
pred <- predict(ridge_opti, s = ridge_lambda, newx = var)
On calcule SCE, SCT et Rcarre_ajuste avec les formules du cours
SCR <- sum((pred - prevision$Temperature_Jour_Suivant)^2)
SCT <- sum((prevision$Temperature_Jour_Suivant - mean(prevision$Temperature_Jour_Suivant))^2)
#on calcule le R² classique
Rcarre <- 1 - (SCR / SCT)
# (SCR)/n-p-1
#R² ajusté = 1 - ____________
# (SCT)/n-p
# SCT = SCR + SCE => SCR = SCT - SCT
RcarreAjust = 1 - ((SCR)/(nrow(prevision) -22 )/ (SCT/(nrow(prevision)-21)))
SCR## [1] 20234.07SCT## [1] 61042.66Rcarre## [1] 0.6685257RcarreAjust## [1] 0.6683228Un R carré de 0.668 est obtenu
## Lasso
lasso <- cv.glmnet(x = var,
y = prevision$Temperature_Jour_Suivant,
alpha = 1,
standardize = TRUE,
nfolds = 5)
lasso_lambda <- lasso$lambda.min
lasso_lambda## [1] 0.04095132predict_lasso <- predict(lasso, s = lasso_lambda, newx = var)
SCR <- sum((predict_lasso - prevision$Temperature_Jour_Suivant)^2)
SCT <- sum((prevision$Temperature_Jour_Suivant - mean(prevision$Temperature_Jour_Suivant))^2)
#on calcule le R² classique
Rcarre <- 1 - (SCR / SCT)
# (SCR)/n-p-1
#R² ajusté = 1 - ____________
# (SCT)/n-p
# SCT = SCR + SCE => SCR = SCT - SCT
RcarreAjust = 1 - ((SCR)/(nrow(prevision) -22 )/ (SCT/(nrow(prevision)-21)))
SCR## [1] 20234.07SCT## [1] 61042.66Rcarre## [1] 0.6685257RcarreAjust## [1] 0.6683228Un R² de 0.668 a été trouvé
Comparaison
Pour les comparaisons qui suivent, nous utilisons nos différents
Test de l’AUC, on utilise un des modèles prédictifs sur le jeu d’entrainement (1324 obs) et sur le jeu de test (331 obs).
# library(ROCR)
# library(pROC)
# #prediction
# PredictEntrainement = predict(reg.mul, newdata = prevision)
# PredictTesting = predict(reg.mul, newdata = testing)
#
# #jeu d'entrainement
# auc(prevision$Temperature_Jour_Suivant, PredictEntrainement)
#
# #jeu de test
# auc(testing$Temperature_Jour_Suivant, PredictTesting)Le résultat est une prédiction parfaite, cependant cet indicateur n’est pas interprétable pour notre modèle.
En effet, l’AUC est valable sur la prédiction de classe , pas sur des prévisions qualitatives.
Comparons nos régressions avec un indicateur ..
Les différents R² ajusté obtenus, sont :
0.99 : Régression sans la constante
0.67 : Régression avec toutes les variables
0.60 : Régression avec variables significatives
0.67 : Régression forward
0.67 : Régression backward
0.67 : Régression sur ACP
0.67 : Régression Ridge
0.67 : Régression Lasso
La régression sans la constante obtiens un R² ajusté de 1, ce qui est logique puisque la droite de régression passe par l’origine.
La régression avec toutes les variables obtient un R² de 0.67, lorsque nous essayons d’obtenir une régression “presque aussi bien” en retirons les variables les moins explicatives, le R² chute à 0.60.
Ensuite, en utilisant les méthodes de protection de la régression, nous obtenons le même R² avec les régressions forward,backward,sur ACP, Ridge que avec la régression avec toutes les variables ! Nous voyons l’efficacité de ces méthodes.
En partant d’une base open-source non construit pour cet effet, nous avons réussis à prédire (à notre échelle) la température du jour suivant en l’exprimant linéairment par rapport à un ensemble de variables explicatives.
## Ouverture (Bonus)
Le test est fait à J+1 : nous essayons de prédire la température du lendemain, mais en changeant notre décalage , nous pouvons prédire à J+n.
Pour un décalage à J+10, on effectue 
Voici les résultats du R² à J+100 par exemple :
J+100
0.99 : Régression sans la constante
0.21 : Régression avec toutes les variables
0.18 : Régression avec variables significatives
0.21 : Régression forward
0.21 : Régression backward
0.21 : Régression sur ACP
0.21 : Régression Ridge
0.21 : Régression Lasso
Les résultats à J+100 sont inexploitable.
Affichage des Rcarres sur 10 jours
require(ggplot2)
require(reshape2)## Le chargement a nécessité le package : reshape2joursDeca <- c( 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 , 10)
Rcarre_Avec_Protection_Regression <- c(0.67 , 0.57 , 0.52 , 0.50 ,0.49 ,0.47 ,0.45 ,0.45 ,0.44 ,0.44)
Rcarre_Avec_Var_Significatives <- c(0.60 , 0.53 , 0.49 , 0.47 ,0.45 ,0.43 ,0.41 ,0.40 ,0.40 ,0.40)
df <- data.frame(joursDeca,Rcarre_Avec_Protection_Regression,Rcarre_Avec_Var_Significatives)
df <- melt(df , id.vars = 'joursDeca', variable.name = 'Rcarres')
ggplot(df, aes(joursDeca,value)) + geom_line(aes(colour = Rcarres))
Les prédictions sont plus difficiles avec le décalage des jours.
Néanmoins, un résultat à J+1 de 0.64 de Rcarre ajusté reste mauvais.
Nous avons fait de nombreuses hypothèses pour expliquer ce résultat.
Nous pensons que les contraintes du jeu de données ont beaucoup influé : les jours où il n’y a aucune observation ne sont pas négligables ainsi que la sélection de la première observation par jour ne rend pas compte de l’ensemble de toutes les données de la journées.
De plus,nous estimons que nos modèle ne tiennent pas en compte la tendance et la saisonnalité (nos données sont chronologique).
Vérifions cette hypothèse en comparant les moyennes de température entre les années.
On utilise pour cela l’Anova
# modele <- aov( TempÃ.rature ~mois ,prevision)
#
# anova(modele)
modele <- aov( Température ~annee ,prevision2)
anova(modele)## Analysis of Variance Table
##
## Response: Température
## Df Sum Sq Mean Sq F value Pr(>F)
## annee 1 227 227.09 6.1726 0.01307 *
## Residuals 1654 60851 36.79
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Test sur les résidus
residus <- residuals(modele)
shapiro.test(residus)##
## Shapiro-Wilk normality test
##
## data: residus
## W = 0.99586, p-value = 0.0001693Le test n’est pas significatif.
Puisque le test n’est pas significatif,Réalisons un test de Kruskal.
# kruskal.test(TempÃ.rature ~mois ,prevision)
# prevision$mois <- NULL
kruskal.test(Température ~annee ,prevision2)##
## Kruskal-Wallis rank sum test
##
## data: Température by annee
## Kruskal-Wallis chi-squared = 51.701, df = 6, p-value = 2.141e-09#un test de bartlett
bartlett.test(Température ~annee ,prevision2)##
## Bartlett test of homogeneity of variances
##
## data: Température by annee
## Bartlett's K-squared = 25.177, df = 6, p-value = 0.0003167La différence est très significative.
La différence de moyenne est aussi significative , même lorsque l’on regroupe les années en 3 groupes.
# #library(questionr)
# #icut(prevision, annee)
#
#
# ## Recodage de prevision$annee en prevision$annee_rec
# prevision$annee_rec <- cut(prevision$annee,
# include.lowest = TRUE,
# right = FALSE,
# dig.lab = 4,
# breaks = 3
# )
Les différences étant très significative, cela explique les résultats peu convaincants à J+1.
Nous ajouterons prochainement un modèle plus adapté à notre problèmatique, tel le modèle de Buys Ballot.