Charger le package mixOmics
ou l’installer si besoin :
library(mixOmics)
# install.packages("mixOmics")
On note Variance expliquée,en \(\%\) de la variance totale de \(X\), dans le cas de la régression PLS
(sparse
et même DA
), la grandeur, pour une matrice de covariables centrée \(X\), un poids associé \(u\)
\[ Var_{X,u}= \frac{||Xu||_2^2}{||X||_F^2}*100, \] où
\(||Xu||_2^2=\{\text{Somme de tous les éléments de $Xu$ élevés au carré}\}\),
\(||X||_F^2=\{\text{Somme de tous les éléments de $X$ élevés au carré}\}\).
Il a été demandé en cours de l’intérêt de construire une composante qui a une faible variance potentielle. L’exemple qui suit montre l’impérieuse nécessité de ne pas se concentrer seulement sur la variance expliquée lorsque l’on fait une analyse en PLS.
On regarde comment l’ajout d’une composante peut avoir un intérêt certain sur les résultats, en analyse discriminante, sur des données simulées.
n1 <- 50
n2 <- 30
n3 <- 20
n4 <- 10
mu <- 10
X1 <- cbind(rnorm(n1,mu),rnorm(n1,0),rnorm(n1,0),rnorm(n1,0))
X2 <- cbind(rnorm(n2,0),rnorm(n2,mu),rnorm(n2,0),rnorm(n2,0))
X3 <- cbind(rnorm(n3,0),rnorm(n3,0),rnorm(n3,mu),rnorm(n3,0))
X4 <- cbind(rnorm(n4,0),rnorm(n4,0),rnorm(n4,mu),rnorm(n4,mu))
X <- scale(rbind(X1,X2,X3,X4))
Y <- as.factor(c(rep("A",n1),rep("B",n2),rep("C",n3),rep("D",n4) ))
On construit un modèle à \(3\) composantes que l’on observe
r <- 3
mm <- plsda(X,Y,ncomp = r)
plot(as.data.frame(mm$variates$X),col=Y,pch=as.numeric(Y)+14)
On peut aussi représenter les densités marginales.
library(ggplot2)
toplot <- data.frame(cbind(mm$variates$X,as.factor(Y)))
alpha <- 0.4
cols = c("black","red","green","blue")
ggplot(toplot, aes(comp.1, fill = Y)) +
geom_density(alpha = alpha)+
scale_fill_manual(values=cols)
ggplot(toplot, aes(comp.2, fill = Y)) +
geom_density(alpha = alpha)+
scale_fill_manual(values=cols)
ggplot(toplot, aes(comp.3, fill = Y)) +
geom_density(alpha = alpha)+
scale_fill_manual(values=cols)
On peut aussi s’attarder à calculer la variance expliquée par chaque axe, en \(\%\) de la variance totale de \(X\):
varTotX <- sum((X)^2)
var_comps <- diag(crossprod(mm$variates$X)/varTotX)*100
barplot(var_comps,pch=16,ylim=c(0,max(var_comps)),ylab="Variance explained (% of total X variance)");abline(h=0)
Ces valeurs sont aussi accessibles via l’argument explained_variance
de tout modèle PLS de mixOmics. On y retrouve aussi les variances calculées pour Y.
Remarque :
Ces variances ne sont pas forcemment décroissantes par rapport à l’indice de la composante construite contrairement à l’ACP. Faire bien attention :)
Question On voit très nettement que les deux premières composantes apportent énormemment d’information. Mais qu’en est-il de la troisième ? Faut-il en conclure que cette composante n’est pas utile ?
errs <- rep(0,r)
for(h in 1:r){
mm <- plsda(X,Y,ncomp = h)
cvcv <- perf(mm,progressBar = F,validation = "loo")
errs[h] <- cvcv$error.rate$overall[h,1]
}
names(errs) <- c("comp 1","comp 2","comp 3")
barplot(errs,pch=16,ylim=c(0,max(errs)),ylab="Error rate");abline(h=0)
Il est clair que ne construire que 2 composante amènerait à faire une grosse erreur potentielle, \(\sim 10\%\) d’erreur en leave-one-out.
La question est maintenant :
Est-ce intéressant de construire cette composante qui permet de réduire drastiquement l’erreur de prédiction mais qui apporte peu de variance ?
La réponse est bien entendu OUI. Se tromper dans \(\sim 10\%\) des cas c’est énorme. D’autant plus que les erreurs sont toutes commises pour un même groupe, le groupe bleu, si l’on ne construisait pas cette troisième composante, légère en terme de variance mais pas inutile du tout.
Ici on voit bien que la variance expliquée seule ne suffit pas pour faire le choix de modèle. La validation croisée s’avère être nécessaire. Ce n’est pas anecdotique, en grande dimension, il faut faire très attention, dans le cas de la PLS, aux conclusions que l’on peut faire en regardant seulement les variances expliquées par chacun des axes. La validation croisée permet seule de répondre à la question initiale de la prédiction du Y.