Modélisation (ou calibration) en grande dimension $p>n$ où les variables explicatives sont les discrétisations de spectres en proche infra rouge (NIR) et la variable à expliquer la part de sucre dans une pâte à gâteau. Différentes méthodes de régression (ridge, lasso, lars, elasticnet, svm, pls, sur composantes principales}...) adaptées à cette situation sont utilisés et leur qualités prédictives comparées. La librairie caret est ensuite utilisée pour industrialiser la stratégie de choix de modèle et de méthode.
Un deuxième scénario sur des données NIR (Tecator) aborde des méthodes non linéaires et l'usage des splines pour lisser et calculer les dérivées.
Le travail présenté s'intéresse à un problème de contrôle de qualité sur une chaîne de fabrication de biscuits (cookies). Il est nécessaire de contrôler le mélange des ingrédients avant cuisson afin de s'assurer que les proportions en lipides, sucre, farine, eau, sont bien respectées c'est-à-dire proches des valeurs nominales de la recette qui a fait la réputation de l'entreprise avant l'industrialisation de la production. Il s'agit de savoir s'il est possible de dépister au plus tôt une dérive afin d'intervenir sur les équipements concernés. Les mesures et analyses, faites dans un laboratoire classique de chimie, sont relativement longues et coûteuses ; elles ne peuvent être entreprises pour un suivi régulier ou même en continue de la production. C'est pour cela qu'il à été décidé l'utilisation d'un spectromètre en proche infrarouge (NIR).
Ce type de problème est très classique en agro-alimentaire et plus généralement dans l'industrie et correspond à une étude dite de calibration en chimiométrie. Une très abondante littérature lui est consacrée.
Les données originales sont dues à Osbone et al. (1984) et ont été souvent utilisées pour la comparaison de méthodes (Stone et al. 1990, Brown et al. 2001, Krämer et al. 2008). Elles sont accessibles dans R au sein du package ppls.
Les mesures ont été faites sur deux échantillons, l'un de taille 40 prévu pour l'apprentissage, l'autre de taille 32 pour les tests. Pour chacun de ces 72 biscuits, les compositions en lipides, sucre, farine, eau, sont mesurées par une approche classique tandis que le spectre est observé sur toutes les longueurs d'ondes entre 1100 et 2498 nanomètres, régulièrement espacés de 2 nanomètres. Nous avons donc 700 valeurs observées, ou variables potentiellement explicatives, par échantillon de pâte à biscuit.
Typiquement, cette étude se déroule dans un contexte de grande dimension avec $p>n$.
L'objectif principal est de répondre à la question suivante : est-il possible de prévoir ces quantités à partir des spectres ? En cas de réponse positive, le gain de temps, donc d'argent, serait immédiat. L'étude est restreinte à la seule modélisation du taux de sucre pour la recherche d'un meilleur modèle de prévision. Il s'agit donc d'évaluer, comparer différentes stratégies et méthodes pour aboutir à celle la plus efficace. Le travail est découpé en deux parties, dans la première sont testées différentes méthodes en construisant les modèles à la main ou pas à pas tandis que la 2ème partie fait appel à une librairie caret qui propose une comparaison automatique d'optimisation des paramètres pour un vaste ensemble de méthodes.
Répondre aux questions Q.
Nicole Krämer, Anne Laure Boulesteix et Gerhard Tutz (2008). Penalized Partial Least Squares with applications to B-spline transformations and functional data, Chemometrics and Intelligent Laboratory Systems 94, 60–69.
Max Kuhn (2008). Building Predictive Models in R Using the caret Package, Journal of Statistical Software, 28, 5.
B. G. Osborne, T. Fearn, A. R. Miller et S. Douglas (1984). Application of Near Infrared Reflectance spectroscopy to the compositional analysis of biscuits and biscuit doughs, J. Sci. Food Agric. 35, 99–105.
In [ ]:
library(ppls)
data(cookie)
# extraire le taux de sucre et les spectres
cook = data.frame(cookie[,702],cookie[,1:700])
names(cook)= c("sucre",paste("X",1:700,sep=""))
In [ ]:
# la variable à expliquer
options(repr.plot.width=4, repr.plot.height=4)
hist(cook[,"sucre"], main="")
Q La figure ci-dessus représente la distribution des valeurs de concentration en sucre. Commentez en quelques mots.
Q Les figures ci-dessous représentent les différents spectres (l'absorbance) en fonction de la fréquence (pas du temps) pour chaque biscuit, puis leur distribution. La couleur est fonction de la concentration de sucre. Que dire de ces courbes, de leur apparence commune, de leur régularité et de l'influence sur les corrélations des variables après discrétisation ?
In [ ]:
# les spectres colorés par l'intensité en sucre
coul=rainbow(20)[as.integer(as.factor(as.integer(cook[,1]-10)))]
ts.plot(t(cook[,-1]),col=coul)
In [ ]:
# Les diagrammes boîtes puis variances des variables explicatives
options(repr.plot.width=4, repr.plot.height=4)
boxplot(cook[,-1])
In [ ]:
hist(var(cook[,-1]), main="")
In [ ]:
acp=prcomp(cook[,-1])
options(repr.plot.width=4, repr.plot.height=3)
plot(acp)
In [ ]:
options(repr.plot.width=5, repr.plot.height=5)
biplot(acp)
In [ ]:
options(repr.plot.width=4, repr.plot.height=4)
plot(acp$x,col=coul,pch=19)
In [ ]:
plot.ts(acp$rotation[,1:10], main="")
Q Quel traitement préalable sous contrôle de l'acp serait envisageable sur ces courbes pour une étude exploratoire? Néanmoins, serait-il pertinent au vu de l'objectif de modélisation de ces données?
Une fonction sera utilisée pour une représentation homogène des résidus des modèles :
In [ ]:
plot.res=function(x,y,titre="titre")
{
plot(x,y,col="blue",ylab="Résidus",xlab="Valeurs predites",main=titre,pch=19)
abline(h=0,col="green")
}
La modélisation est construite méthode par méthode pour mieux en assimiler le déroulement. Chaque méthode de modélisation possède des spécificités, notamment dans la manière d'optimiser plus ou moins facilement les valeurs des paramètres. Il est indispensable de bien se familiariser avec ces différents éléments.
Les 40 premières lignes sont usuellement considérées pour l'apprentissage tandis que les 32 dernières sont l'échantillon test. C'est fait pour éventuellement comparer avec des résultats de la littérature.
Q Compte tenu de l'objectif de modélisation et de la structure des données, peut-on envisager l’estimation d’un modèle linéaire gaussien? Pourquoi?
Q Une approche classique de type "sélection de variables" en régression permettrait-elle de résoudre ce problème? Pourquoi?
In [ ]:
# extraire apprentissage et test
cook.app=cook[1:40,]
cook.test=cook[41:72,]
In [ ]:
library(MASS)
# paramètres en fonction de la pénalisation
plot(lm.ridge(sucre ~ ., data=cook.app,lambda = seq(0,1,0.01)))
cook.ridge=lm.ridge(sucre ~ ., data=cook.app,lambda = seq(0,0.2,0.001))
Q Comment optimiser la valeur de lambda?
Q Quelle stratégie est mise en place ci-dessous? Que fait la fonction choix.kappa?
In [ ]:
# Sélection des indices de validation croisée
library(pls)
set.seed(87)
cvseg=cvsegments(nrow(cook.app),k=4,type="random")
In [ ]:
choix.kappa=function(kappamax,cvseg,nbe=100){
press=rep(0,nbe)
for (i in 1:length(cvseg)){
valid=cook.app[unlist(cvseg[i]),]
modele=lm.ridge(sucre~.,data=cook.app[unlist(cvseg[-i]),],lambda=seq(0,kappamax,length=nbe))
coeff=coef(modele)
prediction=matrix(coeff[,1],nrow(coeff),nrow(valid))+coeff[,-1]%*%t(data.matrix(valid[,-1]))
press=press+rowSums((matrix(valid[,1],nrow(coeff),nrow(valid),byrow=T)-prediction)^2)
}
kappaet=seq(0,kappamax,length=nbe)[which.min(press)]
return(list(kappaet=kappaet,press=press))
}
In [ ]:
# exécution
res=choix.kappa(0.5,cvseg,nbe=1000)
options(repr.plot.width=3, repr.plot.height=3)
plot(seq(0,0.5,length=1000),res$press,pch=19)
Q Comment interpréter le résulat?
In [ ]:
# valeur optimale
kappaet=res$kappaet
cook.ridgeo=lm.ridge(sucre~.,data=cook.app,lambda=kappaet)
coeff=coef(cook.ridgeo)
# Calcul des valeurs ajustées et des résidus
fit.rid=rep(coeff[1],nrow(cook.app))+as.vector(coeff[-1]%*%t(data.matrix(cook.app[,-1])))
plot(fit.rid,cook.app[,"sucre"],pch=19)
Q La régression mettant en oeuvre le paramètre optimal fournit le résultat ci-dessus. Commenter la qualité de l'ajustement. Que vaut approximativement le $R^2$?
In [ ]:
res.rid=fit.rid-cook.app[,"sucre"]
options(repr.plot.width=4, repr.plot.height=4)
plot.res(fit.rid,res.rid,titre="")
Q Commenter le graphique précédent.
Q Quels sont les calculs ci-dessous? Noter le résultat pour comparer avec les autres méthodes.
In [ ]:
ychap=rep(coeff[1],nrow(cook.test))+as.vector(coeff[-1]%*%t(data.matrix(cook.test[,-1])))
options(repr.plot.width=4, repr.plot.height=4)
mean((cook.test[,1]-ychap)^2)
In [ ]:
cook.pls=plsr(sucre~.,ncomp=28,data=cook.app,validation="CV",segments=cvseg)
options(repr.plot.width=4, repr.plot.height=4)
plot(cook.pls,pch=19)
Q Commenter la précédure ci-dessous et le résultat obtenu.
In [ ]:
msepcv.pls=MSEP(cook.pls,estimate=c("train","CV"))
plot(msepcv.pls,type="l")
In [ ]:
ncompo=which.min(msepcv.pls$val["CV",,])-1
# modèle optimal
cook.plso=plsr(sucre~.,ncomp=ncompo,data=cook.app)
# résidus
fit.pls=predict(cook.plso,ncomp=1:ncompo)[,,ncompo]
plot(fit.pls,cook.app[,"sucre"], pch=19)
In [ ]:
res.pls=fit.pls-cook.app[,"sucre"]
plot.res(fit.pls,res.pls)
Q Commenter les deux graphiques ci-dessus.
Q Commenter les résultats ci-dessous.
In [ ]:
ychap=predict(cook.plso,newdata=cook.test)[,1,ncompo]
plot(ychap,cook.test[,1])
mean((cook.test[,"sucre"]-ychap)^2)
Attention, une observation (la 23ème) semble se démarquer de l'échantillon. Ce n'est pas flagrant mais elle est signalée comme atypique (outlier) dans la littérature. La repérer sur le graphique; elle sera prise en compte dans la section suivante.
Trois graphiques analogues à ceux de l'ACP viennent compléter les résultats mais ceux-ci sont moins utiles ou pertinents qu'en version 2 de la PLS quand la variable à expliquer $Y$ est multidimensionnelle.
In [ ]:
scoreplot(cook.pls)
corrplot(cook.pls)
loadingplot(cook.pls)
In [ ]:
cook.pcr=pcr(sucre~.,ncomp=28,data=cook.app,validation="CV",segments=cvseg)
msepcv.pcr=MSEP(cook.pcr,estimate= c("train","CV"))
plot(msepcv.pcr,type="l")
In [ ]:
ncompo=which.min(msepcv.pcr$val["CV",,])-1
cook.pcro=pcr(sucre~.,ncomp=ncompo,data=cook.app)
Q Commenter le résultat ci-dessous.
In [ ]:
fit.pcr=predict(cook.pcro,ncomp=1:6)[,,6]
res.pcr=fit.pcr-cook.app[,"sucre"]
chap=predict(cook.pcro,newdata=cook.test)[,1,ncompo]
mean((cook.test[,"sucre"]-ychap)^2)
In [ ]:
library(randomForest)
# estimation
cook.rf=randomForest(sucre~.,data=cook.app,xtest=cook.test[,-1],ytest=
cook.test[,"sucre"],do.trace=50,mtry=50,corr.bias=TRUE)
pred.rfr=cook.rf$test$predicted
mean((pred.rfr-cook.test[,"sucre"])^2)
In [ ]:
pred.rfr=cook.rf$test$predicted
mean((pred.rfr-cook.test[,"sucre"])^2)
Q Que dire de la qualité de ce modèle non linéaire
In [ ]:
library(lars)
frac.delta=seq(0,1,length=1000)
mse.cv=cv.lars(data.matrix(cook.app[,-1]),cook.app[,1],K=4,se=F,index=frac.delta,use.Gram=F)
In [ ]:
frac.delta.o=frac.delta[which.min(mse.cv$cv)]
cook.lasso=lars(data.matrix(cook.app[,-1]),cook.app[,1],use.Gram=F)
fit.lasso=predict(cook.lasso,data.matrix(cook.app[,-1]),s=frac.delta.o,mode="fraction")$fit
plot(fit.lasso,cook.app[,"sucre"], pch=19)
In [ ]:
pred.lasso=predict(cook.lasso,data.matrix(cook.test[,-1]),s=frac.delta.o,mode="fraction")$fit
mean((pred.lasso-cook.test[,"sucre"])^2)
Q Commenter le résultat précédent.
Q Que représente le graphique ci-dessous?
In [ ]:
plot(cook.lasso,breaks=F)
Kunh (2008) a développé tout un environnement (librairie caret) qui facilite modélisation et optimisation des paramètres pour quasiment toutes les méthodes de classification et/ou régression disponibles dans des librairies R et, il y en a beaucoup. Chacune possède une logique ou une stratégie particulière de gestion des paramètres; cet outil permet d'en uniformiser la syntaxe. Cela concerne la transmission des données, du ou des paramètres de complexité à optimiser lors de l'appel des fonctions et le calcul des prévisions.
La librairie possède plusieurs fonctionnalités qui sont mises en oeuvre et illustrées successivement dans le cas présent de la régression:
Les mêmes fonctionnalités sont proposées en situation de classification supervisée avec des spécificités pour le type de l'erreur (mauvais classement, courbe ROC).
De façon générale, le défi pour la suite, est d'arriver à faire mieux en terme de qualité de prévision que la régression PLS.
Osbone et al. (1984) et leur successeurs s'accordent sur le fait que l'observation 23 est atypique. Ils font remarquer de plus que les extrémités des spectres sont bruitées ; en les retirant, les erreurs de prévision sont réduites. Enfin, Goutils et Fearn (1996) suggèrent également que la discrétisation est trop fine. Que devient la prévision avec quatre fois moins de variables ?
Il s'agit donc de reprendre les calculs des erreurs de prévisions avec l'extraction ci-dessous qui supprime la 23ème observation, les 50 premières et dernières valeurs du spectre tout en ne considérant qu'une valeur sur 4:
In [ ]:
cook2=cook[-23,c(1,seq(50,650,4))]
Q Que font les commandes suivantes?
In [ ]:
inTrain = 1:39
trainDescr=cook2[inTrain,-1]
testDescr=cook2[-inTrain,-1]
trainY=cook2[inTrain,1]
testY=cook2[-inTrain,1]
Q Intérêt des commandes suivantes?
In [ ]:
library(doParallel)
cl <- makeCluster(4)
registerDoParallel(cl)
Q Que calculent les transformations suivantes? Purquoi est-ce important pour une méthode comme la régression PLS?
In [ ]:
library(caret)
xTrans=preProcess(trainDescr)
trainDescr=predict(xTrans,trainDescr)
testDescr=predict(xTrans,testDescr)
Q Que définissent les commandes suivantes?
In [ ]:
btControl=trainControl(method="boot",number=100)
cvControl=trainControl(method="cv",number=4)
In [ ]:
# régression ridge
set.seed(2)
ridgeFit = train(trainDescr, trainY,method = "ridge", tuneLength = 8,trControl = cvControl)
ridgeFit
In [ ]:
options(repr.plot.width=4, repr.plot.height=3)
plot(ridgeFit)
In [ ]:
# régression pls
set.seed(2)
plsFit = train(trainDescr, trainY,method = "pls", tuneLength = 12,trControl = cvControl)
plsFit
In [ ]:
plot(plsFit)
In [ ]:
# régression lasso
set.seed(2)
lassoFit = train(trainDescr, trainY,method = "lasso", tuneLength = 8,trControl = cvControl)
lassoFit
In [ ]:
plot(lassoFit)
In [ ]:
# régression lars
set.seed(2)
larsFit = train(trainDescr, trainY,method = "lars", tuneLength = 8,trControl = cvControl)
larsFit
In [ ]:
plot(larsFit)
Q Quelles est plus précisément la méthode ci-dessous? Que sont ses paramètres?
In [ ]:
# elasticnet
set.seed(2)
enetFit = train(trainDescr, trainY,method = "enet", tuneLength = 8,trControl = cvControl)
enetFit
In [ ]:
plot(enetFit)
In [ ]:
# support vector machine
set.seed(2)
svmFit = train(trainDescr, trainY,method = "svmLinear",trControl = cvControl)
svmFit
p
In [ ]:
models=list(ridge=ridgeFit,pls=plsFit,lasso=lassoFit,lars=larsFit,elasticnet=enetFit,svm=svmFit)
testPred=predict(models, newdata = testDescr)
lapply(testPred, function(x) mean((x-testY)^2))
resPlot=extractPrediction(models, testX=testDescr, testY=testY)
plotObsVsPred(resPlot)
In [ ]:
# Attention: exécuter la fonction en annexe
models=c("ridge","pls","lasso","lars","enet","svmLinear")
noptim=c(10,10,10,10,5,5)
# Initialiser le générateur et fixer le nombre d’itérations
# Changer ces valeurs. Attention au temps de calcul! Être patient!
Niter=3 ; Init=11
# Appel de la fonction définie en annexe
X=cook2[,-1]; Y= cook2[,1]
pred.cook=pred.autom(X,Y,methodes=models,N=Niter,xinit=Init,size=noptim)
In [ ]:
# Calcul des risques
prev=pred.cook$pred
obs=pred.cook$obs
# carrés des différences
dif=lapply(prev,function(x)(x-obs)^2)
# risque pour chaque échantillon
moy=lapply(dif,function(x)apply(x,2,mean))
options(repr.plot.width=6, repr.plot.height=4)
boxplot(data.frame(moy))
Q Quel est le graphique ci-dessus. Comment l'interpréter?
In [ ]:
lapply(moy,mean)
In [ ]:
pred.autom=function(X,Y,p=1/2,methodes=c("knn","rf"),size=c(10,2),xinit=11,N=10,typerr="cv",number=4,type="raw") {
# Fonction de prévision de N échantillons tests
# par une liste de méthodes de régression
# ou classification (uniquement 2 classes)
# Optimisation des paramètres par validation
# croisée (défaut) ou bootstrap ou... (cf. caret)
# X : matrice ou frame des variables explicatives
# Y : variable cible quantitative ou qualitative
# p : proportion entre apprentissage et test
# methodes : liste des méthodes de rdiscrimination
# size : e grille des paramètres à optimiser
# xinit : générateur de nombres aléatoires
# N : nombre de réplications apprentissage/test
# typerr : "cv" ou "boo" ou "oob"
# number : nombre de répétitions CV ou bootstrap
# pred : liste des matrices de prévision
# type d’erreur
Control=trainControl(method=typerr,number=number)
# initialisation du générateur
set.seed(xinit)
# liste de matrices stockant les prévisions
# une par méthode
inTrain=createDataPartition(Y,p=p,list=FALSE)
ntest=length(Y[-inTrain])
pred=vector("list",length(methodes))
names(pred)=methodes
pred=lapply(pred,function(x)x=matrix(0,nrow=ntest,ncol=N))
obs=matrix(0,ntest,N)
set.seed(xinit)
for(i in 1:N) {
# N itérations
# indices de l’échantillon d’apprentissage
inTrain=createDataPartition(Y,p=p,list=FALSE)
# Extraction des échantillons
trainDescr=X[inTrain,]
testDescr=X[-inTrain,]
trainY=Y[inTrain]
testY=Y[-inTrain]
# stockage des observés de testY
obs[,i]=testY
# centrage et réduction des variables
xTrans=preProcess(trainDescr)
trainDescr=predict(xTrans,trainDescr)
testDescr=predict(xTrans,testDescr)
# estimation et optimisation des modèles
# pour chaque méthode de la liste
for(j in 1:length(methodes)) {
# modélisation
modFit = train(trainDescr, trainY,method = methodes[j], tuneLength = size[j],trControl = Control)
# prévisions
if (type=="prob") pred[[j]][,i]=predict(modFit,newdata = testDescr,type=type)[,1]
else pred[[j]][,i]=predict(modFit,
newdata = testDescr)
}}
list(pred=pred,obs=obs)
# résultats
}