1 TP ozone : Modèle linéaire gaussien, binomial, courbe ROC

pour répondre au problème de prévision. Ceci passe par la mise en place d’un
TP ozone : Modèle linéaire gaussien, protocole très strict afin de s’assurer d’un minimum d’objectivité pour cette
comparaison.
binomial, courbe ROC
Penser à conserver dans des fichiers les commandes successives utilisées
ainsi que les principaux résultats tableaux et graphes.
Résumé Toutes les opérations sont réalisées dans R avec l’appui de bibliothèques
complémentaires éventuellement à télécharger (mlbench, MASS, boot,
Comparaison sur le même jeu de données des qualités de prévision
class, e1071, rpart, nnet, ipred, gbm, randomForest).
de plusieurs modèles obtenus par analyse de la covariance et régres-
D’autres logiciels sont bien évidemment utilisables (Splus, SAS, SPSS...).
sion logistique. Utilisation d’un échantillon test et courbe ROC.
Néanmoins, dans le cas de SAS, cela nécessiterait l’accès au module Enter-
prise Miner pour certaines méthodes (arbres, neurones) encore peu diffusé
1 Introduction car très cher. De plus certaines procédures du module classique SAS/Stat ne
comportent pas d’option de sélection automatique de variables quantitatives et
1.1 Objectif qualitatives.
L’objectif, sur ces données, est d’améliorer la prévision déterministe
1.2 Rappel : protocole de comparaison
(MOCAGE), calculée par les services de MétéoFrance, de la concentration
d’ozone dans certaines stations de prélèvement. Il s’agit d’un problème dit La démarche mise en œuvre enchaîne les étapes suivantes :
d’adaptation statistique d’une prévision locale de modèles à trop grande 1. Après une éventuelle première étape descriptive uni ou multidimension-
échelle en s’aidant d’autres variables également prévues par MétéoFrance, nelle visant à repérer les incohérences, les variables non significatives
mais à plus petite échelle (température, force du vent...). Plus précisé- ou de distribution exotique, les individus non concernés. . . et à étudier les
ment, deux variables peuvent être prévues : soit la concentration quantitative structures des données, procéder à un tirage aléatoire d’un échantillon test
d’ozone, soit le dépassement (qualitatif) d’un certain seuil fixé à 150µg 1 . Dans qui ne sera utilisé que lors de la dernière étape.
chaque cas, deux approches sont considérées : soit prévoir la concentration
quantitative puis en déduire l’éventuel dépassement ou bien prévoir directe- 2. Sur la partie restante qui sera découpée en échantillon d’apprentissage
ment le dépassement. Dans le premier cas, il s’agit d’abord d’une régression (des paramètres du modèle) et échantillon de validation (pour estimation
tandis que dans le deuxième il s’agit d’un problème de discrimination ou de sans biais du taux de mauvais classés), optimiser les choix afférents à
régression logistique. Question : quelles sont les meilleures méthodes et stra- chacune des méthodes de régression et/ou discrimination :
tégies pour prévoir la concentration d’ozone du lendemain d’une part et l’oc- • variables et interactions à prendre en compte dans la régression linéaire
curence d’un pic de pollution d’autre part ? ou logistique,
• variables et méthode pour l’analyse discriminante,
On se propose de tester différentes méthodes : régression logistique, ana- • nombre de nœuds dans l’arbre de régression ou de classification,
lyse discriminante, réseau de neurones, arbre de décision, agrégation d’arbres • architecture (nombre de couches, de neurones par couche, fonctions de
(bagging, boosting, random forest), SVM. L’objectif final, à ne pas perdre de transferts, nombre de cycles. . . ) du perceptron,
vue, est la comparaison de ces méthodes afin de déterminer la plus efficace • algorithme d’agrégation,
• noyau et complexité des SVMs.
1. Le seuil légal a été porté à 180µg mais celui-ci est bien trop rarement dépassé par les
observations et serait trop difficile à prévoir par estimation d’un modèle.
 2 TP ozone : Modèle linéaire gaussien, binomial, courbe ROC

Remarques : # Vérification du contenu

• En cas d’échantillon petit face au nombre des variables il est recom- summary(ozone)
mandé d’itérer la procédure de découpage par validation croisée, c’est # Changement du type de la variable jour
le cas de ces données, afin de réduire la variance des estimations des ozone[,"JOUR"]=as.factor(ozone[,"JOUR"])
erreurs de classement.
Remarquer le type des variables. Il est nécessaire d’en étudier la distribution.
3. Comparaison finale des qualités de prévision sur la base du taux de mal
classés pour le seul échantillon test qui est resté à l’écart de tout effort ou hist(ozone[,"O3obs"]);hist(ozone[,"MOCAGE"])
"acharnement" pour l’optimisation des modèles. hist(ozone[,"TEMPE"]);hist(ozone[,"RMH2O"])
• Attention, ne pas "tricher" en modifiant le modèle obtenu lors de l’étape hist(ozone[,"NO2"]);hist(ozone[,"NO"])
précédente afin d’améliorer le résultat sur l’échantillon test ! hist(ozone[,"VentMOD"]);hist(ozone[,"VentANG"])
• Le critère utilisé dépend du problème : erreur quadratique, taux de mau-
Des transformations sont proposées pour rendre certaines distributions plus
vais classement, AUC (aire sous la courbe ROC)...
symétriques et ainsi plus "gaussiennes".
1.3 Prise en charge des données ozone[,"SRMH2O"]=sqrt(ozone[,"RMH2O"])
Les données ont été extraites et mises en forme par le service concerné de ozone[,"LNO2"]=log(ozone[,"NO2"])
MétéoFrance. Elles sont disponibles dans le fichier ozone.dat et décrites ozone[,"LNO"]=log(ozone[,"NO"])
par les variables suivantes :
Vérifier l’opportunité de ces transformations puis retirer les variables initiales
JOUR Le type de jour ; férié (1) ou pas (0) ; et construire la variable "dépassement de seuil".
O3obs La concentration d’ozone effectivement observée le lendemain à 17h ozone=ozone[,c(1:4,8:13)]
locales correspondant souvent au maximum de pollution observée ; ozone[,"DepSeuil"]=as.factor(ozone[,"O3obs"]>150)
MOCAGE Prévision de cette pollution obtenue par un modèle déterministe
de mécanique des fluides (équation de Navier et Stockes) ; pour obtenir le fichier qui sera effectivement utilisé.
TEMPE Température prévue par MétéoFrance pour le lendemain 17h ; Il est également intéressant de regarder les corrélations entre les différentes
variables.
RMH2O Rapport d’humidité ;
NO2 Concentration en dioxyde d’azote ; cor(ozone[,-c(1,5,11)])
pairs(ozone)
NO Concentration en monoxyde d’azote ;
STATION Lieu de l’observation : Aix-en-Provence, Rambouillet, Munchhau- 1.4 Extraction des échantillons apprentissage et test
sen, Cadarache et Plan de Cuques ;
VentMOD Force du vent ; Le programme ci-dessous réalise l’extraction du sous-ensemble des données
d’apprentissage et de test. Attention, chaque étudiant ou binôme tire un échan-
VentANG Orientation du vent. tillon différent ; il est donc "normal" de ne pas obtenir les mêmes modèles.
# Lecture des données Utiliser trois chiffres au hasard, en remplacement de "XXX" ci-dessous,
ozone=read.table("ozone.dat",header=T) comme initialisation du générateur de nombres aléatoires.
 3 TP ozone : Modèle linéaire gaussien, binomial, courbe ROC

set.seed(XXX) # initialisation du générateur 2.1 Modèle linéaire

# Extraction des échantillons
test.ratio=.2 # part de l’échantillon test
Le modèle de régression linéaire simple intégrant des variables qualitatives
npop=nrow(ozone) # nombre de lignes dans les données correspondant est estimé avec la fonction aov mieux adaptée à l’analyse de
nvar=ncol(ozone) # nombre de colonnes covariance.
# taille de l’échantillon test
ntest=ceiling(npop*test.ratio) # estimation du modèle sans interaction
# indices de l’échantillon test reg.lm=aov(O3obs~.,data=datappr)
testi=sample(1:npop,ntest) summary(reg.lm)
# indices de l’échantillon d’apprentissage # Extraction des résidus et des valeurs ajustées
appri=setdiff(1:npop,testi) res.lm=reg.lm$residuals
fit.lm=reg.lm$fitted.values
Construction des échantillons pour la régression (prévision de la concentra- # Graphe des résidus
tion). plot(fit.lm,res.lm)
# Définition de 2 fonctions pour un graphe coloré
# construction de l’échantillon d’apprentissage # et des échelles fixes sur les axes
datappr=ozone[appri,-11] plot.res=function(x,y,titre="titre")
# construction de l’échantillon test {
datestr=ozone[testi,-11] plot(x,y,col="blue",xlim=c(0,300),ylim=c(-100,100),
summary(datappr) # vérifications ylab="Résidus",xlab="Valeurs predites",main=titre)
summary(datestr) # points(x2,y,col="red")
abline(h=0,col="green")
Construction des échantillons pour la discrimination (dépassement du seuil). }
# construction de l’échantillon d’apprentissage plot.res(fit.lm,res.lm,"")
datappq=ozone[appri,-2] plot.fit=function(x,y,titre="titre"){
# construction de l’échantillon test plot(x,y,col="blue",xlim=c(0,400),ylim=c(0,400),
datestq=ozone[testi,-2] ylab="Ajustées",xlab="Observations O3",main=titre)
summary(datappq) # vérifications abline(0,1,col="green")
summary(datestq) abline(150,0,col="pink")
abline(v=150,col="pink")}
plot.fit(datappr[,"O3obs"],fit.lm,
2 Prévision par modèle gaussien "Ajustees modele lineaire")
# Graphe des résidus au modèle MOCAGE
Le premier modèle à tester est un simple modèle de régression linéaire mais, plot.res(datappr[,"MOCAGE"],
comme certaines variables sont qualitatives, il s’agit d’une analyse de cova- datappr[,"MOCAGE"]-datappr[,"O3obs"],"")
riance. D’autre part, on s’intéresse à savoir si des interactions sont à prendre # Regroupement des graphiques sur une même page
en compte. Le modèle devient alors polynomiale d’ordre 2 ou quadratique. par(mfrow=c(2,2))
 4 TP ozone : Modèle linéaire gaussien, binomial, courbe ROC

hist(res.lm) 2.3 Sélection de modèle par sélection de variables

qqnorm(res.lm)
qqline(res.lm,col=2) Afin de faire cette sélection de variables, la library MASS est utilisée ainsi
plot.res(fit.lm,res.lm) que les critères d’Akaike et BIC. Les méthodes de sélection ascendante, des-
plot.fit(datappr[,"O3obs"],fit.lm, cendante et stepwise sont mise en oeuvre.
"Ajustees modele lineaire") 2.3.1 Sélection par AIC et backward
Contrôler les graphes des résidus, que conclure de ce modèle ? L’étude suivante library(MASS)
met en œuvre toutes les interactions d’ordre 2 pour ajuster le modèle. modlinAIC_B=stepAIC(reg.lm,~.,trace=TRUE,
On peut s’intéresser à la qualité d’apprentissage de ce modèle linéaire com- direction=c("backward"))
plet en calculant son erreur quadratique moyenne d’apprentissage. summary(modlinAIC_B)
mean(res.lm^2) Noter que certains paramètres sont significatifs selon la méthode de sélec-
tion descendante d’AIC mais ne le sont pas selon le test de Fisher.
2.2 Nouvelle paramétrisation 2.3.2 Sélection par AIC et forward
Afin de faciliter l’interprétation des résultats concernant les variables quali- mod0=lm(O3obs~1,data=datappr)
tatives, on peut introduire une nouvelle paramétrisation à l’aide de contrastes. modlinAIC_F=stepAIC(mod0,O3obs~JOUR+MOCAGE+TEMPE+
Par défaut, la référence est prise par rapport à la première modalité de la va- STATION+VentMOD+VentANG+SRMH2O+LNO2+LNO,trace=TRUE,
riable qualitative. Les paramètres indiqués pour les variables Jour et Station direction=c("forward"))
indiquent l’écart estimé par rapport à cette référence. Il est plus itéressant de summary(modlinAIC_F)
se référer à la moyenne des observations sur toutes les modalités des variables
qualitatives et d’interpréter les coefficients comme des écarts à cette moyenne. 2.3.3 Sélection par AIC et stepwise
contrasts(datappr$JOUR)=contr.sum( modlinAIC_S=stepAIC(reg.lm,~.,trace=TRUE,
levels(datappr$JOUR)) direction=c("both"))
contrasts(datappr$STATION)=contr.sum( #L’option "both" est l’option par défaut
levels(datappr$STATION)) #de la commande StepAIC
summary(modlinAIC_S)
modlin=lm(O3obs~.,datappr)
summary(modlin) 2.3.4 Sélection par BIC et stepwise
# Attention aux noms des variables : JOUR1 = 0
On peut effectuer cette sélection de variable à l’aide du critère de BIC en
# JOUR2 =1, STATION1 = Aix , etc...
remplaçant la pénalité 2 par le log du nombre de données.
# La somme des coefficients associée aux variables
qualitatives est nulle. Ainsi les coefficients des #k=log(npop) pour BIC au lieu de k=2 pour AIC
modalités JOUR=1 et STATION=Ramb ne sont pas modlinBIC_S=stepAIC(reg.lm,~.,trace=TRUE,
affichés, car ils sont l’opposé de la somme des direction=c("both"),k=log(npop))
paramètres estimés des autres modalités. summary(modlinBIC_S)
 5 TP ozone : Modèle linéaire gaussien, binomial, courbe ROC

Puisque la pénalité est plus forte, on constate qu’effectivement le modèle sé- # du critère BIC par méthode descendante
lectionné par BIC est plus parcimonieux que celui sélectionné par AIC. reg.glm.stepBIC_D=step(reg.glm,
direction="backward",k=log(npop))
2.4 Calcul de l’erreur d’apprentissage anova(reg.glm.stepBIC_D,test="F")
Avant de vérifier la qualité de prédiciton de chacun des modèles construits, # Extraction des valeurs ajustées et des résidus
il faut vérifier leur qualité d’apprentissage. Les différents graphes sont donc à fit.glmBIC_D=reg.glm.stepBIC_D$fitted.values
tracer et les erreurs quadratiques moyennes d’apprentissage sont à comparer. res.glmBIC_D=reg.glm.stepBIC_D$residuals
# Graphe des résidus
#Modèle stepwise AIC plot.res(fit.glmBIC_D,res.glmBIC_D,"Residus
mean((predict(modlinAIC_S)-datappr$O3obs)^2) modele quadratique BIC_D")

#Modèle stepwise BIC On remarque que la présence de certains interactions ou variables sont perti-
mean((predict(modlinBIC_S)-datappr$O3obs)^2) nentes au sens du critère d’Akaïke mais pas significatives au sens du test de Fi-
sher. Cette présence dans le modèle peut être plus finement analysée en consi-
2.5 Modèle quadratique dérant une estimation de l’erreur par validation croisée. L’idée est de retirer
une à une les variables ou interactions les moins significatives et de voir com-
Le modèle de régression quadratique est estimé avec la fonction glm qui ment se comporte la validation croisée. D’autre part, si la procédure pas-à-pas
permet aussi une sélection automatique de modèle. Les méthodes descendante, conduit à un modèle différent, l’estimation de l’erreur par validation croisée
ascendante et stepwise des critères AIC et BIC sont à mettre en oeuvre. Afin de permet également d’optimiser le choix.
ne pas faire trop de répétitions dans l’énoncé de ce TP, la méthode descendante
pour AIC et BIC sont présentées à l’aide de la commande step du package stats. L’estimation des erreurs par validation croisée est calculée en utilisant une
fonction proposée dans la bibliothèque boot.
# Estimation du modèle avec interactions d’ordre 2
reg.glm=glm(O3obs~(.)^2,data=datappr) library(boot) # chargement de la bibliothèque
# validation croisée 10-plis
# Recherche du meilleur modèle au sens # modèle complet
# du critère d’Akaïke par méthode descendante cv.glm(datappr, reg.glm, K=10)$delta[1]
reg.glm.stepAIC_D=step(reg.glm, # modèle "Akaïke" méthode descendante
direction="backward") cv.glm(datappr, reg.glm.stepAIC_D, K=10)$delta[1]
anova(reg.glm.stepAIC_D,test="F")
La première commande suivante permet de récupérer le modèle optimal avant
# Extraction des valeurs ajustées et des résidus
de redéfinir celui-ci en retirant l’interaction la moins significative avant de ré-
fit.glmAIC_D=reg.glm.stepAIC_D$fitted.values
estimer l’erreur par validation croisée dans le cas du modèle sélectionné par
res.glmAIC_D=reg.glm.stepAIC_D$residuals
Akaike méthode descendante.
# Graphe des résidus
plot.res(fit.glmAIC_D,res.glmAIC_D,"Residus #la liste des effets du modèle optimal
modele quadratique AIC_D") reg.glm.stepAIC_D$formula
cv.glm(datappr, glm(O3obs ~JOUR + MOCAGE + TEMPE +
# Recherche du meilleur modèle au sens STATION + VentMOD + SRMH2O + LNO2 + LNO+
 6 TP ozone : Modèle linéaire gaussien, binomial, courbe ROC

JOUR:LNO + MOCAGE:VentMOD + MOCAGE:LNO2 + F=table[2,1]/(table[2,1]+table[1,1])

MOCAGE:LNO + TEMPE:STATION + TEMPE:SRMH2O + PSS=H-F
TEMPE:LNO2 + TEMPE:LNO + STATION:VentMOD + print(PSS)
STATION:SRMH2O + STATION:LNO2 + }
STATION:LNO + VentMOD:SRMH2O,
data=datappr),K=10)$delta[1] 2.6.2 Comparaison des qualités prédictives sur l’échantillon test
On souhaite retenir le meilleur modèle obtenu pour prédire l’échantillon test
2.6 Prévision de l’échantillon test et estimer ainsi sans biais une erreur de prévision. Deux erreurs sont estimées ;
2.6.1 Quelques scores de prévision la première est celle de l’erreur quadratique moyenne dans le cas de la régres-
sion tandis que la deuxième est issue de la matrice de confusion qui croise
On note A l’occurrence du phénomène d’intérêt et N A sa non-occurrence. les dépassements de seuils prédits avec ceux effectivement observés. Les cal-
culs sont présentés pour le modèle linéaitre avec interaction sélection par la
Le taux de succès global est défini méthode d’Akaike descendante.
Observé par
a+d # Calcul des prévisions
Prévu A N A .
A a b a+b+c+d pred.glmAIC_D=predict(reg.glm.stepAIC_D,
newdata=datestr)
NA c d Ce n’est pas forcément un bon in-
# Erreur quadratique moyenne de prévision
dicateur si A est rare.
mean((pred.glmAIC_D-datestr[,"O3obs"])^2)
Soient H le taux de bonnes prévisions et F le taux de fausses alertes. Ces deux # Matrice de confusion pour la prévision du
taux se calculent de la manière suivante # dépassement de seuil
table(pred.glmAIC_D>150,datestr[,"O3obs"]>150)
a b pss(datestr[,"O3obs"]>150,pred.glmAIC_D>150)
H= et F = .
a+c b+d
Noter ces erreurs pour les comparer avec celles obtenues par les autres mé-
Grâce à ces deux taux, on peut déterminer le score de Pierce thodes. Par exemple avec la prévision de MOCAGE.
P SS = H − F ∈ [−1; 1]. # matrice de confusion du modèle MOCAGE
table(datestq[,"MOCAGE"]>150,datestq[,"DepSeuil"])
Le cas P SS = 0 correspond à une prévision climatologique : prévision au pss(datestq[,"DepSeuil"],datestq[,"MOCAGE"]>150)
hasard qui respecte les fréquences des événements observés.
Si P SS > 0, alors le modèle apporte de l’information. On peut créer une
fonction permettant de calculer le score de Pierce. 3 Prévision par modèle binomial
pss=function(obs,prev) Plutôt que de prévoir la concentration puis le dépassement, on peut se poser
{ la question de savoir s’il ne serait pas pertinent de prévoir directement la pré-
table=table(prev,obs) sence ou l’absence d’un dépassement. La variable à modéliser étant binaire,
H=table[2,2]/(table[2,2]+table[1,2]) c’est la régression logistique qui va être employée. Comme pour la régression,
 7 TP ozone : Modèle linéaire gaussien, binomial, courbe ROC

différentes stratégies de choix de modèle peuvent être utilisées et comparées

avant d’estimer l’erreur de prévision sur l’échantillon test. # algorithme backward
log.qm=glm(DepSeuil~(.)^2,data=datappq,
3.1 Régression logistique sans interaction family=binomial)
log.qm.step2=step(log.qm,direction="backward",
# estimation du modèle complet
family=binomial)
log.lm=glm(DepSeuil~.,data=datappq,family=binomial) # significativité des paramètres
# significativité des paramètres
anova(log.qm.step2,test="Chisq")
anova(log.lm,test="Chisq")
# Recherche d’un modèle optimal au sens d’Akaïke Deux modèles sont en concurrence, il s’agit de les comparer en minimisant
log.lm.step=step(log.lm,direction="backward") l’erreur calculée par validation croisée. La même procédure est utilisée que
anova(log.lm.step,test="Chisq") pour le modèle linéaire classique mais une autre fonction de coût adaptée à une
# matrice de confusion de l’échantillon variable binaire doit être précisée. Elle est définie ci-dessous puis appliquée au
# d’apprentissage et erreur apparente calcul de l’erreur apparente de prévision.
table(log.lm.step$fitted.values>0.5,
datappq[,"DepSeuil"]) cout = function(r, pi=0)
mean(abs(as.integer(r)-pi)>0.5)
# Application de cette fonction pour
3.2 Régression logistique avec interaction
# l’erreur apparente ou d’ajustement :
# estimation du modèle complet cout(as.integer(datappq$DepSeuil)-1,
log.qm=glm(DepSeuil~(.)^2,data=datappq, log.qm.step1$fitted.values)
family=binomial) cout(as.integer(datappq$DepSeuil)-1,
log.qm.step2$fitted.values)
Avec autant de paramètres, l’algorithme de régression peut rencontrer quelques
soucis si la régression ajuste "trop bien" les données. On propose alors une Estimation de l’erreur de prévision par validation croisée.
optimisation du modèle par une méthode pas-à-pas à comparer avec le modèle library(boot)
fournit par l’algorithme descendant. cv.glm(datappq, log.qm.step1, cout,K=10)$delta[1]
cv.glm(datappq, log.qm.step2, cout,K=10)$delta[1]
# régression avec le modèle minimum
log.qm=glm(DepSeuil~1,data=datappq,family=binomial) Retenir le "meilleur" modèle.
# algorithme stepwise en précisant le plus grand
# modèle possible 3.3 Prévision de l’échantillon test
log.qm.step1=step(log.qm,direction="both", 3.3.1 Score de Brier
scope=list(lower=~1,upper=~(JOUR + MOCAGE +
TEMPE + STATION + VentMOD + VentANG + LNO2 + En plus du score de Pierce, dans le cas d’un modèle binomial, on peut égale-
LNO + SRMH2O)^2), family=binomial) ment s’intéresser au score de Brier. Ce score évalue les prévisions probabilistes
# significativité des paramètres de variables ne pouvant prendre qu’un nombre fini de valeurs (ou "catégories").
anova(log.qm.step1,test="Chisq") C’est le cas de la variable DepSeuil. Si on note pi la prévision de la probabilité
 8 TP ozone : Modèle linéaire gaussien, binomial, courbe ROC

de réalisation de l’événement Dépassement de seuil de la i-ème observation library(verification)

et oi l’observation de la modalité de l’observation i (valant 1 si l’évènement DepSeuil=as.numeric(datestq[,"DepSeuil"])-1
Dépassement de seuil est réalisé et 0 sinon), alors le score de Brier, noté BS, brier(DepSeuil,pred.log)
se calcule de la façon suivante sur un échantillon de taille n :
Mémoriser les résultats obtenus pour comparer avec les autres méthodes.
n
1X 2 3.3.3 Diagramme de fiabilité
BS = (pi − oi ) .
n i=1
Les commandes suivantes permettent de réaliser le diagramme de fiabi-
lité d’un modèle prédisant une probabilité de réalisation d’un événement. Ce
Cette formule peut se généraliser au cas d’une variable cible qualitative à
graphe représente les probabilités prédites par le modèle en fontion de la fré-
m ≥ 2 modalités de la façon suivante :
quence d’observation de l’événement en question.
n m
1 XX 2 plot.fiabilite=function(n,pred,obs,mod){
BS = (pi,j − oi,j ) ,
2n i=1 j=1
fiab=rep(0,n)
où oi,j vaut 1 si la i-ème observation est de la catégorie j et 0 sinon ; et prob=rep(0,n)
où pi,j désigne la propabilité prévue pour la j-ème catégorie pour la i-ème for (i in 1:n){
observation. fiab[i]=mean(obs[(pred<=(i/n))
& (pred>((i-1)/n))])
prob[i]=mean(pred[(pred<=(i/n))
Le score de Brier est une distance dans le domaine des probabilités. Aussi, & (pred>((i-1)/n))])
la valeur de ce score sera d’autant plus faible que la prévision sera bonne et }
une prévision parfaite obtiendra un score de 0. Le plus mauvais score est de 1.
La commande plot(prob,fiab,xlab="Probabilités prédites",
xlim=c(0,1),ylim=c(0,1),
brier ylab="Fréquences",
main=mod,pch=8,col="red",cex=1.3)
du package verification, peut être utilisée.
abline(a=0,b=1,col="orange")
3.3.2 Calcul sur l’échantillon test }

pred.log=predict(log.qm.step1,newdata=datestq, DepSeuilObs=datestq[,"DepSeuil"]
type="response") DepSeuilObsF=as.integer(DepSeuilObs)-1
# Matrice de confusion pour la prévision du plot.fiabilite(10,pred.log,DepSeuilObsF,
# dépassement de seuil "Diagramme de fiabilité - Régression Logistique")
table(pred.log>0.5,datestq[,"DepSeuil"])
La commande suivante fournit exactement le même graphique mais avec des
#Calcul des scores de Pierce et de Brier commandes plus synthétiques.
pss(datestq[,"DepSeuil"],pred.log>0.5) A<- verify(DepSeuilObs, pred.log,
 9 TP ozone : Modèle linéaire gaussien, binomial, courbe ROC

frcst.type = "prob", obs.type = "binary") perflogit=performance(predlogit, "tpr","fpr")

reliability.plot(A, plot(perflogit,col=1)
main="Diagramme de fiabilité - predlog.AIC_D") aire1=performance(predlogit,"auc")
aire1@y.values
Notre modèle aura une bonne qualité de prédiction si les points représentés
sur le diagramme de fiabilité sont proches de la première bissectrice. Il est également possible de construire une courbe ROC en association de la
prévision obtenue à partir d’un modèle gaussien. En effet, la variation du seuil
Afin de juger quelles sont les erreurs commises par notre modèle, il est pos- théorique de dépassement (150) va faire varier les proportions respectives des
sible, sur un même graphique, de tracer les probabilités prédites par le modèle taux de vrais et faux positifs. Cela revient encore à faire varier le seuil d’une
et la modalité observée, ie dans notre cas la variable dépassement de seuil. "proba" pour les valeurs de prévisions divisées par 300.

plot.prob=function(pred,obs,mod){ rocglm=pred.glmAIC_D/300
predglm=prediction(rocglm,datestq[,"DepSeuil"])
rg=sort(pred,index.return=TRUE) perfglm=performance(predglm, "tpr","fpr")
plot(rg$x,type="l",col="orange",lwd=3, plot(perfglm,col=2,add=TRUE)
xlab="Observations",ylab="Probabilités", aire2=performance(predglm,"auc")
main=mod) aire2@y.values
points(obs[rg$ix],col="red",cex=0.3,pch=1)
Les résultats obtenus dépendent évidemment en plus de l’échantillonnage
legend(x="topleft",legend=c("Probabilités
initial entre apprentissage et test. Dans le cas où les courbes se croisent, cela
prédites","Dépassement de seuil de
signifie qu’il n’y a pas de prévision uniformément meilleure de l’occurrence
pollution"),lty=c(1,0),
de dépassement. Cela dépend de la sensibilité ou de la spécificité retenue pour
pch=c(-1,1),col=c("orange","red"),
le modèle. Ceci souligne l’importance de la bonne définition du critère à utili-
cex=0.8)
ser pour le choix d’une "meilleure" méthode. Ce choix dépend directement de
abline(h=0.5)
celui , "politique" ou "économique" de sensibilité et / ou spécificité du modèle
abline(v=150)
retenu. En d’autres termes, quel taux de fausse alerte, avec des imputations
}
économiques évidentes, est supportable au regard des dépassements non dé-
plot.prob(pred.log,DepSeuilObs,
tectés et donc de la dégradation sanitaire de la population à risque ?
"Régression Logistique")
C’est une fois ce choix arrêté que le statisticien peut opérer une comparaison
3.3.4 Courbe ROC et aire AUC des méthodes en présence.

Pour finir, la courbe ROC et l’aire AUC peuvent être utilisés afin de compa-
rer l’ensemble des modèles construits.
library(ROCR)
roclogit=predict(log.qm.step2,newdata=datestq,
type="response")
predlogit=prediction(roclogit,datestq[,"DepSeuil"])