Modélisation (ou calibration) en grande dimension. Les variables explicatives sont les discrétisations de spectres en proche infra rouge (NIR) et la variable à expliquer la part de matière grasse dans des échantillons de viande. Différente méthodes de régression lasso, krls, pls, svm, mars...) sont utilisées et leur qualités prédictives comparées. La prise en compte de la forme des spectres par l'introduction des dérivées après lissage spline est utile sur ce jeu de données. La librairie caret est directement utilisée pour industrialiser la stratégie de choix de modèle et de méthode. Ce scénario vient compléter celui des données NIR de pâte à biscuit.
These data are recorded on a Tecator Infratec Food and Feed Analyzer working in the wavelength range 850 - 1050 nm by the Near Infrared Transmission (NIT) principle. Each sample contains finely chopped pure meat with different moisture, fat and protein contents.}\footnote{The data are available in the public domain with no responsability from the original data source. The data can be redistributed as long as this permission note is attached.
Chaque échantillon de viande est caractérisé par des mesures d'absorption pour 100 longueurs d'ondes échantillonnées dans le proche infra rouge ainsi que par des mesures chimiques de la proportion (pourcentage) en eau, graisse et protéine. Les mesures d'absorption ont déjà été transformées ($\log_{10}$). Nous nous attacherons à modéliser le pourcentage de matière grasse.
Ces données ont été largement étudiées dans la littérature car elles servent de jeu de donnée test (benchmark) pour comparer des prétraitements des données (centrage, réduction sur données brutes ou après ACP, pls, décomposition sur base de splines et dérivations, réduction de dimension par ACP ou PLS) des méthodes de modélisation (régression linéaire, svm, réseaux de neurones, cart...), des techniques de sélection de variables (stepwise, par critère bayésien, informaiton mutuelle, par algorithme stochastique...) et toutes les combinaisons de ces différentes approches...
Contrairement à la plupart des données NIR du domaine public, ces données présentent une composante non linéaire certes légère mais présente. De plus, la prise en compte des dérivées ou tout du moins des différences d'ordre 2 des spectres est aussi pertinente, d'où l'intérêt qu'elles ont sucité et suscitent toujours dans la littérature du métier, mais aussi celle statistique ou de machine learning, pour justifier l'utilisation de nouvelles méthodes: non-linéaires et/ou fonctionnelles.
Les données sont accessibles dans R au sein du package caret. Elles sont contenues dans deux matrices: absorp de $n=215$ lignes et $p=100$ colonnes, les absorptions par longueur d'onde, et endpoints de $n=215$ lignes et $p=3$ colonnes, les pourcentages d'eau, graisses, et protéines. Les 129 premières observations sont généralement utilisées comme échantillon d'apprentissage.
L'objectif principal est la recherche d'un meilleur modèle de prévision du pourcentage de matière grasse. Il s'agit d'évaluer, comparer, différentes stratégies et méthodes pour aboutir à celle la plus efficace. Le travail est découpé en trois parties. La première explore les spectres, leur lissage et de leur dérivées; la deuxième modélise les spectres initiaux tandis que la suivante prend en compte les dérivées des spectre après lissage spline avant finalement d'automatiser les traitements pour optimiser les choix.
Trois variables très coorélées entre elles peuvent être modélisées et prédites à partir des mesutes d'absorbtion ou spectres.
In [ ]:
library(caret)
# Variables à expliquer
data(tecator)
summary(endpoints)
In [ ]:
# Spectre des absorptions ou variables explicatives
dim(absorp)
In [ ]:
# corrélation des variables à expliquer
options(repr.plot.width=4, repr.plot.height=4)
splom(endpoints)
In [ ]:
## tracer un échantillon aléatoire de 10 spectres
# sélection des observations
set.seed(1)
inSubset = sample(1:dim(endpoints)[1], 10)
absorpSubset = absorp[inSubset,]
endpointSubset = endpoints[inSubset, 2]
# tri des spectres selon la première valeure
newOrder = order(absorpSubset[,1])
absorpSubset = absorpSubset[newOrder,]
endpointSubset = endpointSubset[newOrder]
# définition des couleurs
plotColors = rainbow(10)
# tracé des échelles
plot(absorpSubset[1,], type = "n", ylim = range(absorpSubset), xlim = c(0, 105), xlab = "Wavelength Index", ylab = "Absorption")
# tracé des spectres et du taux de graisse
for(i in 1:10){
points(absorpSubset[i,],type = "l",
col = plotColors[i], lwd = 2)
text(105,absorpSubset[i,100], endpointSubset[i],
col = plotColors[i])
}
title("Profiles for 10 Random Samples")
La moyenne de chaque spectre ne semble pas être une information pertinente pour l'estimation du taux de graisse. En revanche, ils présentent des formes remarquable dans certaines zones, d'où l'intérêt éventuel d'étudier la dérivée de ces courbes.
In [ ]:
# la variable à expliquer choisie, taux de matières grasses, est la deuxième
gras=endpoints[,2]
hist(gras)
In [ ]:
# les diagrammes boîtes des variables explicatives
options(repr.plot.width=8, repr.plot.height=4)
boxplot(data.frame(absorp))
In [ ]:
acp=prcomp(absorp)
options(repr.plot.width=3, repr.plot.height=3)
plot(acp)
In [ ]:
options(repr.plot.width=4, repr.plot.height=4)
biplot(acp)
In [ ]:
plot(acp$x)
In [ ]:
plot.ts(acp$rotation[,1:10])
Très grande importance de la première composante principale. Elle a pour forme celle moyenne des courbes dont les décalages génèrent une variance importante.
Choix du paramètre de lissage ou degré de liberté qui est la trace de la hat matrix ou matrice de lissage. Le premier spectre, son lissage et ses dérivées sont tracés pour différentes valeurs du paramètre {\tt df}. Si celui-ci n'est pas fourni (dernier cas), il est optimisé par validation croisée afin de reconstruire au mieux la courbe.
In [ ]:
library(splines)
options(repr.plot.width=6, repr.plot.height=6)
par(mfcol=c(3,3))
degre=c(3,10,20)
for (i in 1:3){
model=smooth.spline(x=1:100,y=absorp[1,],df=degre[i])
lisse=predict(model)$y
deriv1=predict(model,deriv=1)$y
deriv2=predict(model,deriv=2)$y
ts.plot(cbind(lisse,absorp[1,]))
ts.plot(deriv1)
ts.plot(deriv2)
}
In [ ]:
model=smooth.spline(x=1:100,y=absorp[1,])
lisse=predict(model)$y
deriv1=predict(model,deriv=1)$y
deriv2=predict(model,deriv=2)$y
options(repr.plot.width=6, repr.plot.height=3)
par(mfcol=c(1,3))
ts.plot(cbind(lisse,absorp[1,]))
ts.plot(deriv1)
ts.plot(deriv2)
Il semblerait raisonnable de lisser un peu plus que ce que suggère la validation croisée afin d'obtenir des dérivées plus lisses. L'idéal serait évidemment d'optimiser ce paramètre en cherchant à prévoir au mieux le pourcentage de matières grasses.
Génération des bases de B-splines cubiques et tracer de quelques unes d'entre elles.
In [ ]:
base=bs(1:100,df=model$fit$nk)
options(repr.plot.width=3, repr.plot.height=3)
ts.plot(base[,c(1,20,30,64)])
Pour chaque spectre, les commandes ci-dessous calculent un lissage des spectres par B-splines cubiques : noeuds équirépartis, degré fixé à 15. Le programme fournit les dérivées première et seconde aux points de discrétisation, les coefficients dans la base de base de B-Splines.
In [ ]:
degre=15
# initialisation des matrices
coeff=matrix(ncol=model$fit$nk,nrow=215)
lisse=matrix(ncol=100,nrow=215)
deriv1=matrix(ncol=100,nrow=215)
deriv2=matrix(ncol=100,nrow=215)
# itération sur tous les spectres
for (obs in 1:215){
model=smooth.spline(x=1:100,y=absorp[obs,],df=degre)
coeff[obs,]=model$fit$coef
lisse[obs,]=predict(model)$y
deriv1[obs,]=predict(model,deriv=1)$y
deriv2[obs,]=predict(model,deriv=2)$y
}
In [ ]:
dimnames(absorp)[[2]]=1:100
In [ ]:
# Indices usuelles des échantillons
# d'apprentissage et de test
inTrain = 1:129
trainDescr=absorp[inTrain,]
testDescr=absorp[-inTrain,]
trainY=gras[inTrain]
testY=gras[-inTrain]
Il est recommandé de centrer et réduire les variables dans plusieurs méthodes. C'est fait systématiquement et simplement en utilisant évidemment les mêmes transformations sur l'échantillon test que celles mises en place sur l'échantillon d'apprentissage. D'autres part, l'erreur de prévision est estimée par validation croisée.
In [ ]:
xTrans=preProcess(trainDescr)
trainDescr=predict(xTrans,trainDescr)
testDescr=predict(xTrans,testDescr)
cvControl=trainControl(method="cv",number=10)
In [ ]:
library(doParallel)
cl <- makeCluster(4)
registerDoParallel(cl)
In [ ]:
#1 Régression linéaire leapBackward
set.seed(2)
llmFit = train(trainDescr, trainY,method = "leapBackward", tuneLength = 20,trControl = cvControl)
plot(llmFit)
In [ ]:
#2 régression pls
set.seed(2)
plsFit = train(trainDescr, trainY,method = "pls", tuneLength = 30,trControl = cvControl)
plot(plsFit)
In [ ]:
#3 Ridge regression
set.seed(2)
ridgeFit = train(trainDescr, trainY,method = "ridge", tuneLength = 8,trControl = cvControl)
plot(ridgeFit)
In [ ]:
#4 Regression elastic net
set.seed(2)
enetFit = train(trainDescr, trainY,method = "enet", tuneLength = 8,trControl = cvControl)
plot(enetFit)
In [ ]:
#5 Support vector machine noyau linéaire
set.seed(2)
svmFit = train(trainDescr, trainY, method = "svmLinear",trControl = cvControl)
In [ ]:
#6 Support vector machine noyau gaussien
set.seed(2)
svmgFit = train(trainDescr, trainY,method = "svmRadial", tuneLength = 8,trControl = cvControl)
plot(svmgFit)
In [ ]:
#8 Kernel-based Regularized Least Squares (KRLS)
set.seed(2)
krlsFit = train(trainDescr, trainY,method = "krlsRadial", tuneLength = 10,trControl = cvControl)
plot(krlsFit)
In [ ]:
models=list(leaplm=llmFit,pls=plsFit,elasticnet=enetFit,svm=svmFit,svmg=svmgFit,krls=krlsFit)
testPred=predict(models, newdata = testDescr)
lapply(testPred,function(x)sqrt(mean((x-testY)^2)))
resPlot=extractPrediction(models, testX=testDescr, testY=testY)
plotObsVsPred(resPlot)
Remarquer quels sont les modèle prenant en compte la non linéarité des données, ceux acceptant des valeurs tests relativements atypiques.
In [ ]:
# Mêmes indices de l'échantillon d'apprentissage
dimnames(deriv1)[[2]]=1:100
inTrain = 1:129
trainDescr1=deriv1[inTrain, ]
testDescr1 =deriv1[-inTrain,]
# pré-process
xTrans1=preProcess(trainDescr1)
trainDescr1=predict(xTrans1,trainDescr1)
testDescr1 =predict(xTrans1,testDescr1 )
Le calcul est restreint au 3 meilleures méthodes trouvées précédemment.
In [ ]:
# Régression PLS
set.seed(2)
plsFit1 = train(trainDescr1, trainY,method = "pls", tuneLength = 30,trControl = cvControl)
plot(plsFit1)
In [ ]:
# Kernel-based Regularized Least Squares (KRLS)
set.seed(2)
krlsFit1 = train(trainDescr1, trainY,method = "krlsRadial", tuneLength = 10,trControl = cvControl)
In [ ]:
plot(krlsFit1)
In [ ]:
# Caclul des erreurs et tracé des graphes
models=list(pls1=plsFit1,krls1=krlsFit1)
testPred1=predict(models, newdata = testDescr1)
lapply(testPred1,function(x)sqrt(mean((x-testY)^2)))
resPlot=extractPrediction(models, testX=testDescr1, testY=testY)
plotObsVsPred(resPlot)
L'échantillon est de faible taille, les estimations des erreurs très dépendantes de l'échantillon test sont sujettes à caution et on peut s'interroger sur la réalité des différences entre les différentes méthodes. En regardant les graphiques, on observe qu'il suffit d'une observation pour influencer les résultats. Il est donc important d'itérer le processus sur plusieurs échantillons tests afin de décrire les distributions des erreurs.
Evidemment le temps de calcul s'en ressent.
Exécuter la fonction en annexe avant d'exécuter:
In [ ]:
models=c("pls","krlsRadial")
noptim=c(10,10,10,10)
# Initialiser le générateur et fixer le nombre d’itérations
# Changer ces valeurs. Attention au temps de calcul! Être patient!
Niter=30 ; Init=11
# Appel de la fonction définie en annexe
X=deriv1; Y= gras
pred.tecat=pred.autom(X,Y,methodes=models,N=Niter,xinit=Init,size=noptim)
Puis calculer :
In [ ]:
# Calcul des risques
prev=pred.tecat$pred
obs=pred.tecat$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))
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
}
In [ ]: