ozone <- read.table("../donnees/ozone_complet.txt", header = T, sep = ";")
dim(ozone)[1] 1464 2310 Comparaison des différentes méthodes, étude de cas réels
Préliminaires
Importation de l’ozone
Élimination des individus avec une valeur manquante
indNA <- which(is.na(ozone), arr.ind = T)[,1]
ozone2 <- ozone[-indNA,]Importation du fichier d’ozone sans valeurs manquantes avec les projections
ozone <- read.table("../donnees/ozone_transf.txt", header = T, sep = ";")et préparation du data-frame qui contiendra les résultats de chaque méthode
RES <- data.frame(Y = ozone$maxO3)Pour le moment il ne contient qu’une seule colonne avec les données à prévoir.
Méthodes et comparaison
Nous séparons en 10 blocs en 2 étapes: création des affectations des individus à chaque bloc
nbbloc <- 10
blocseq <- rep(1:nbbloc, length = nrow(ozone))Puis nous utilisons une permutation aléatoire de ces affectations
set.seed(1234)
bloc <- sample(blocseq)
set.seed(1234)
bloc <- sample(blocseq)Régression multiple
Chargement du package pour la sélection de variables
library(leaps)Evaluation de la qualité prédictive de la régression linéaire et de la sélection de variables via BIC (algorithme exhaustif)
for(i in 1:nbbloc){
###MCO global
reg <- lm(maxO3~.,data=ozone[bloc!=i,])
RES[bloc==i,"MCO"] <- predict(reg,ozone[bloc==i,])
###MCO choix
recherche <- regsubsets(maxO3~., int=T, nbest=1, nvmax=22,
data=ozone[bloc!=i,])
resume <- summary(recherche)
nomselec <- colnames(resume$which)[
resume$which[which.min(resume$bic),] ][-1]
formule <- formula(paste("maxO3~",paste(nomselec,collapse="+")))
regbic <- lm(formule,data=ozone[bloc!=i,])
RES[bloc==i,"choix"] <- predict(regbic,ozone[bloc==i,])
}Lasso, ridge et elasticnet
Chargement du package pour lasso, ridge et elasticnet et création des matrices nécessaires à son utilisation :
library(glmnet)
ozone.X <- model.matrix(maxO3~.,data=ozone)[,-1]
ozone.Y <- ozone[,"maxO3"]Évaluation de la qualité prédictive des régressions lasso, ridge et elasticnet:
for(i in 1:nbbloc){
XA <- ozone.X[bloc!=i,]
YA <- ozone.Y[bloc!=i]
XT <- ozone.X[bloc==i,]
###ridge
tmp <- cv.glmnet(XA,YA,alpha=0)
mod <- glmnet(XA,YA,alpha=0,lambda=tmp$lambda.min)
RES[bloc==i,"ridge"] <- predict(mod,XT)
###lasso
tmp <- cv.glmnet(XA,YA,alpha=1)
mod <- glmnet(XA,YA,alpha=0,lambda=tmp$lambda.min)
RES[bloc==i,"lasso"] <- predict(mod,XT)
###elastic
tmp <- cv.glmnet(XA,YA,alpha=0.5)
mod <- glmnet(XA,YA,alpha=.5,lambda=tmp$lambda.min)
RES[bloc==i,"elastic"] <- predict(mod,XT)
}Régressions sur composantes
Chargement du package pour les régressions sur composantes
library(pls)Évaluation de la qualité prédictive des régressions PCR et PLS
for(i in 1:nbbloc){
#####PLS
tmp <- plsr(maxO3~.,data=ozone[bloc!=i,],ncomp=20,
validation="CV",scale=TRUE)
mse <- MSEP(tmp,estimate=c("train","CV"))
npls <- which.min(mse$val["CV",,])-1
mod <- plsr(maxO3~.,ncomp=npls,data=ozone[bloc!=i,],scale=TRUE)
RES[bloc==i,"PLS"] <- predict(mod,ozone[bloc==i,],ncomp=npls)
#####PCR
tmp <- pcr(maxO3~.,data=ozone[bloc!=i,],ncomp=20,
validation="CV",scale=TRUE)
mse <- MSEP(tmp,estimate=c("train","CV"))
npcr <- which.min(mse$val["CV",,])-1
mod <- pcr(maxO3~.,ncomp=npcr,data=ozone[bloc!=i,],scale=TRUE)
RES[bloc==i,"PCR"] <- predict(mod,ozone[bloc==i,],ncomp=npcr)
}Les résultats :
RES |>
dplyr::summarise(across(-Y,~mean((Y-.)^2))) MCO choix ridge lasso elastic PLS PCR
1 187.2726 188.8491 187.8304 187.1023 187.0389 187.2622 187.2685Pour aller plus loin
Création d’un data-frame pour les résultats :
res_rep <- data.frame(MCO=0,choix=0,ridge=0,lasso=0,elastic=0)Régression linéaire
La fonction
sse_reg <- function(don,bloc,b) {
m_reg <- lm(maxO3~.,data=don[bloc!=b,])
previsions <- predict(m_reg,don[bloc==b,])
return(sum((don[bloc==b,"maxO3"]-previsions)^2))
}La qualité de la modélisation
set.seed(1234)
ssereg <- rep(0,20)
for (r in 1:20) {
bloc <- sample(blocseq)
for(b in 1:nbbloc){
ssereg[r] <- ssereg[r] + sse_reg(ozone,bloc,b)
}
}
res_rep$MCO <- round(mean(ssereg/nrow(ozone)),2)Choix de variables
La fonction
library(leaps)
sse_regbic <- function(don,bloc,b,nvmax,method) {
recherche <- regsubsets(maxO3~., int=T, nbest=1,data=don[bloc!=b,],
nvmax=nvmax,method=method)
resume <- summary(recherche)
nomselec <- colnames(resume$which)[resume$which[which.min(resume$bic),]][-1]
formule <- formula(paste("maxO3 ~", paste(nomselec, collapse = "+")))
m_reg <- lm(formule,data=don[bloc!=b,])
previsions <- predict(m_reg,don[bloc==b,])
return(sum((don[bloc==b,"maxO3"]-previsions)^2))
}La qualité de la modélisation
set.seed(1234)
sseregbic <- rep(0,20)
for (r in 1:20) {
bloc <- sample(blocseq)
for(b in 1:nbbloc){
sseregbic[r] <- sseregbic[r] + sse_regbic(ozone,bloc,b,22,"exhaustive")
}
}
res_rep$choix <- mean(sseregbic/nrow(ozone))Lasso
La fonction
library(glmnet)
sse_glmnet <- function(X,Y,bloc,b,a) {
rech <- cv.glmnet(X[bloc!=b,], Y[bloc!=b,drop=FALSE], alpha=a)
prev <- predict(rech, newx=X[bloc==b,], s=rech$lambda.min)
return(sum((Y[bloc==b,"maxO3"] - as.vector(prev))^2))
}La qualité de la modélisation
X <- model.matrix(maxO3~.,data=ozone)[,-1]
Y <- data.matrix(ozone[,"maxO3",drop=FALSE])
set.seed(1234)
sselasso <- rep(0,20)
for (r in 1:20) {
bloc <- sample(blocseq)
for(b in 1:nbbloc){
sselasso[r] <- sselasso[r] + sse_glmnet(X,Y,bloc,b,a=1)
}
}
res_rep$lasso <- round(mean(sselasso/nrow(ozone)),2)RIDGE
La qualité de la modélisation
X <- model.matrix(maxO3~.,data=ozone)[,-1]
Y <- data.matrix(ozone[,"maxO3",drop=FALSE])
set.seed(1234)
sseridge <- rep(0,20)
for (r in 1:20) {
bloc <- sample(blocseq)
for(b in 1:nbbloc){
sseridge[r] <- sseridge[r] + sse_glmnet(X,Y,bloc,b,a=0)
}
}
res_rep$ridge <- mean(sseridge/nrow(ozone))Elastic-net
La qualité de la modélisation
X <- model.matrix(maxO3~.,data=ozone)[,-1]
Y <- data.matrix(ozone[,"maxO3",drop=FALSE])
set.seed(1234)
sseelasticnet <- rep(0,20)
for (r in 1:20) {
bloc <- sample(blocseq)
for(b in 1:nbbloc){
sseelasticnet[r] <- sseelasticnet[r] + sse_glmnet(X,Y,bloc,b,a=0.5)
}
}
res_rep$elastic <- mean(sseelasticnet/nrow(ozone))Les résultats
res_rep MCO choix ridge lasso elastic
1 188.36 189.6025 188.6624 187.92 187.8242