1. Importer les données

don <- read.table("https://r-stat-sc-donnees.github.io/spe_bretagne.txt", sep=";", header=TRUE, row.names=1)
moy.ligne <- apply(don[,-1], 1, mean)
don[,-1] <- sweep(don[,-1], 1, moy.ligne, FUN="-")
et.ligne <- apply(don[,-1], 1, sd)
don[,-1] <- sweep(don[,-1], 1, et.ligne, FUN="/")

2. Représenter les données

hist(don[,"CO"], prob=TRUE, main="", xlab="Teneur en carbone organique")
lines(density(don[,"CO"]))
rug(don[,"CO"], col="red")

which(don[,1]>8)
[1] 79
max(don[-79,1])
[1] 4.89
don <- don[-79,]
coul <- as.numeric(cut(don[,1], quantile(don[,1], prob=seq(0,1,by=1/7)),
include.lowest = TRUE))
palette(terrain.colors(7,alpha=0.5))
lity <- coul
matplot(x=400:2500, y=t(as.matrix(don[,-1])), type="l", lty=1,
col=coul, xlab="Longueur d’onde", ylab="Réflectance")

3. Effectuer une régression PLS après avoir choisi le nombre de composantes PLS

library(pls)

Attachement du package : 㤼㸱pls㤼㸲

The following object is masked from 㤼㸱package:stats㤼㸲:

    loadings
modele.pls <- plsr(CO~., ncomp=100, data=don,scale=TRUE,validation="CV")
msepcv.pls <- MSEP(modele.pls, estimate=c("train","CV"))
palette("default") ## retour aux couleurs habituelles
plot(msepcv.pls, lty=1, type="l", legendpos="topright", main="")

ncomp.pls <- which.min(msepcv.pls$val["CV",,])-1
ncomp.pls
4 comps 
      4 
reg.pls <- plsr(CO~., ncomp=ncomp.pls, data=don, scale=TRUE)

4. Analyser les résidus

res.pls <- residuals(reg.pls)
plot(res.pls[,,ncomp.pls],pch=15,cex=.5,ylab="Résidus",main="")
abline(h=c(-2,0,2), lty=c(2,1,2))

5. Prévoir une nouvelle valeur

donN <- read.table("https://r-stat-sc-donnees.github.io/spe_nouveau.txt", sep=";", header=TRUE, row.names=1)
moy.ligneN <- apply(donN[,-1], 1, mean)
donN[,-1] <- sweep(donN[,-1], 1, moy.ligneN, FUN="-")
318
[1] 318
et.ligneN <- apply(donN[,-1], 1, sd)
donN[,-1] <- sweep(donN[,-1], 1, et.ligneN, FUN="/")
pred <- predict(reg.pls, ncomp=ncomp.pls, newdata=donN[,-1])
donN[,1]
[1] 1.08 1.60 1.85

Pour aller plus loin

plot(reg.pls, plottype="coef", comps=1:ncomp.pls, main="",
legendpos="topleft", xlab="Longueur d’onde", lty=1,labels=400:2500)

colo <- rep(1,nrow(don))
colo[substr(rownames(don),0,4)=="rmqs"] <- 2
plot(reg.pls, plottype="scores", comps=c(1,2), col=colo, asp=1)
abline(h=0,lty=2)
abline(v=0,lty=2)

coul <- gray(seq(0,.9,len=ncol(don)))
plot(modele.pls,plottype="correlation",comps=1:2,col=coul,pch=20)

plot(modele.pls,plottype="loadings", labels=400:2500, comps=1:ncomp.pls,
legendpos="topright",lty=1, xlab="Longueur d’onde", ylab="Loadings")
abline(h=0,lty=2)

LS0tDQp0aXRsZTogIlJlZ3Jlc3Npb24gUExTIg0KYXV0aG9yOiAiSHVzc29uIGV0IGFsLiINCmRhdGU6ICIwNS8wOS8yMDE4Ig0Kb3V0cHV0Og0KICBodG1sX25vdGVib29rOg0KICAgIHRvYzogeWVzDQogICAgdG9jX2RlcHRoOiAzDQogICAgdG9jX2Zsb2F0OiB5ZXMNCiAgaHRtbF9kb2N1bWVudDoNCiAgICB0b2M6IHllcw0KICAgIHRvY19kZXB0aDogJzMnDQogICAgdG9jX2Zsb2F0OiB5ZXMNCi0tLQ0KIyAxLiBJbXBvcnRlciBsZXMgZG9ubsOpZXMNCg0KYGBge3J9DQpkb24gPC0gcmVhZC50YWJsZSgiaHR0cHM6Ly9yLXN0YXQtc2MtZG9ubmVlcy5naXRodWIuaW8vc3BlX2JyZXRhZ25lLnR4dCIsIHNlcD0iOyIsIGhlYWRlcj1UUlVFLCByb3cubmFtZXM9MSkNCmBgYA0KDQpgYGB7cn0NCm1veS5saWduZSA8LSBhcHBseShkb25bLC0xXSwgMSwgbWVhbikNCmRvblssLTFdIDwtIHN3ZWVwKGRvblssLTFdLCAxLCBtb3kubGlnbmUsIEZVTj0iLSIpDQpldC5saWduZSA8LSBhcHBseShkb25bLC0xXSwgMSwgc2QpDQpkb25bLC0xXSA8LSBzd2VlcChkb25bLC0xXSwgMSwgZXQubGlnbmUsIEZVTj0iLyIpDQpgYGANCg0KIyAyLiBSZXByw6lzZW50ZXIgbGVzIGRvbm7DqWVzDQpgYGB7cn0NCmhpc3QoZG9uWywiQ08iXSwgcHJvYj1UUlVFLCBtYWluPSIiLCB4bGFiPSJUZW5ldXIgZW4gY2FyYm9uZSBvcmdhbmlxdWUiKQ0KbGluZXMoZGVuc2l0eShkb25bLCJDTyJdKSkNCnJ1Zyhkb25bLCJDTyJdLCBjb2w9InJlZCIpDQpgYGANCg0KYGBge3J9DQp3aGljaChkb25bLDFdPjgpDQptYXgoZG9uWy03OSwxXSkNCmBgYA0KDQpgYGB7cn0NCmRvbiA8LSBkb25bLTc5LF0NCmBgYA0KDQpgYGB7cn0NCmNvdWwgPC0gYXMubnVtZXJpYyhjdXQoZG9uWywxXSwgcXVhbnRpbGUoZG9uWywxXSwgcHJvYj1zZXEoMCwxLGJ5PTEvNykpLA0KaW5jbHVkZS5sb3dlc3QgPSBUUlVFKSkNCnBhbGV0dGUodGVycmFpbi5jb2xvcnMoNyxhbHBoYT0wLjUpKQ0KbGl0eSA8LSBjb3VsDQptYXRwbG90KHg9NDAwOjI1MDAsIHk9dChhcy5tYXRyaXgoZG9uWywtMV0pKSwgdHlwZT0ibCIsIGx0eT0xLA0KY29sPWNvdWwsIHhsYWI9Ikxvbmd1ZXVyIGTigJlvbmRlIiwgeWxhYj0iUsOpZmxlY3RhbmNlIikNCmBgYA0KDQojIDMuIEVmZmVjdHVlciB1bmUgcsOpZ3Jlc3Npb24gUExTIGFwcsOocyBhdm9pciBjaG9pc2kgbGUgbm9tYnJlIGRlIGNvbXBvc2FudGVzIFBMUw0KDQpgYGB7cn0NCmxpYnJhcnkocGxzKQ0KbW9kZWxlLnBscyA8LSBwbHNyKENPfi4sIG5jb21wPTEwMCwgZGF0YT1kb24sc2NhbGU9VFJVRSx2YWxpZGF0aW9uPSJDViIpDQpgYGANCg0KYGBge3J9DQptc2VwY3YucGxzIDwtIE1TRVAobW9kZWxlLnBscywgZXN0aW1hdGU9YygidHJhaW4iLCJDViIpKQ0KcGFsZXR0ZSgiZGVmYXVsdCIpICMjIHJldG91ciBhdXggY291bGV1cnMgaGFiaXR1ZWxsZXMNCnBsb3QobXNlcGN2LnBscywgbHR5PTEsIHR5cGU9ImwiLCBsZWdlbmRwb3M9InRvcHJpZ2h0IiwgbWFpbj0iIikNCmBgYA0KYGBge3J9DQpuY29tcC5wbHMgPC0gd2hpY2gubWluKG1zZXBjdi5wbHMkdmFsWyJDViIsLF0pLTENCm5jb21wLnBscw0KYGBgDQoNCmBgYHtyfQ0KcmVnLnBscyA8LSBwbHNyKENPfi4sIG5jb21wPW5jb21wLnBscywgZGF0YT1kb24sIHNjYWxlPVRSVUUpDQpgYGANCg0KIyA0LiBBbmFseXNlciBsZXMgcsOpc2lkdXMNCg0KYGBge3J9DQpyZXMucGxzIDwtIHJlc2lkdWFscyhyZWcucGxzKQ0KcGxvdChyZXMucGxzWywsbmNvbXAucGxzXSxwY2g9MTUsY2V4PS41LHlsYWI9IlLDqXNpZHVzIixtYWluPSIiKQ0KYWJsaW5lKGg9YygtMiwwLDIpLCBsdHk9YygyLDEsMikpDQpgYGANCg0KIyA1LiBQcsOpdm9pciB1bmUgbm91dmVsbGUgdmFsZXVyDQoNCmBgYHtyfQ0KZG9uTiA8LSByZWFkLnRhYmxlKCJodHRwczovL3Itc3RhdC1zYy1kb25uZWVzLmdpdGh1Yi5pby9zcGVfbm91dmVhdS50eHQiLCBzZXA9IjsiLCBoZWFkZXI9VFJVRSwgcm93Lm5hbWVzPTEpDQptb3kubGlnbmVOIDwtIGFwcGx5KGRvbk5bLC0xXSwgMSwgbWVhbikNCmRvbk5bLC0xXSA8LSBzd2VlcChkb25OWywtMV0sIDEsIG1veS5saWduZU4sIEZVTj0iLSIpDQozMTgNCmV0LmxpZ25lTiA8LSBhcHBseShkb25OWywtMV0sIDEsIHNkKQ0KZG9uTlssLTFdIDwtIHN3ZWVwKGRvbk5bLC0xXSwgMSwgZXQubGlnbmVOLCBGVU49Ii8iKQ0KYGBgDQoNCmBgYHtyfQ0KcHJlZCA8LSBwcmVkaWN0KHJlZy5wbHMsIG5jb21wPW5jb21wLnBscywgbmV3ZGF0YT1kb25OWywtMV0pDQpgYGANCg0KYGBge3J9DQpkb25OWywxXQ0KYGBgDQoNCiMgUG91ciBhbGxlciBwbHVzIGxvaW4NCmBgYHtyfQ0KcGxvdChyZWcucGxzLCBwbG90dHlwZT0iY29lZiIsIGNvbXBzPTE6bmNvbXAucGxzLCBtYWluPSIiLA0KbGVnZW5kcG9zPSJ0b3BsZWZ0IiwgeGxhYj0iTG9uZ3VldXIgZOKAmW9uZGUiLCBsdHk9MSxsYWJlbHM9NDAwOjI1MDApDQpgYGANCg0KYGBge3J9DQpjb2xvIDwtIHJlcCgxLG5yb3coZG9uKSkNCmNvbG9bc3Vic3RyKHJvd25hbWVzKGRvbiksMCw0KT09InJtcXMiXSA8LSAyDQpwbG90KHJlZy5wbHMsIHBsb3R0eXBlPSJzY29yZXMiLCBjb21wcz1jKDEsMiksIGNvbD1jb2xvLCBhc3A9MSkNCmFibGluZShoPTAsbHR5PTIpDQphYmxpbmUodj0wLGx0eT0yKQ0KYGBgDQoNCmBgYHtyfQ0KY291bCA8LSBncmF5KHNlcSgwLC45LGxlbj1uY29sKGRvbikpKQ0KcGxvdChtb2RlbGUucGxzLHBsb3R0eXBlPSJjb3JyZWxhdGlvbiIsY29tcHM9MToyLGNvbD1jb3VsLHBjaD0yMCkNCmBgYA0KDQpgYGB7cn0NCnBsb3QobW9kZWxlLnBscyxwbG90dHlwZT0ibG9hZGluZ3MiLCBsYWJlbHM9NDAwOjI1MDAsIGNvbXBzPTE6bmNvbXAucGxzLA0KbGVnZW5kcG9zPSJ0b3ByaWdodCIsbHR5PTEsIHhsYWI9Ikxvbmd1ZXVyIGTigJlvbmRlIiwgeWxhYj0iTG9hZGluZ3MiKQ0KYWJsaW5lKGg9MCxsdHk9MikNCmBgYA0KDQo=