######################################################### # # Correction Tp 5 # ########################################################## rm(list=ls()) # -------------------------------------------------------- # 1. Recuperation des donnees # -------------------------------------------------------- # On sauvegarde le jeu olympic.txt au bon endroit. # Ce sont les perfomances de 34 athletes aux 10 epreuves du decathlon. # http://pbil.univ-lyon1.fr/ade4/ade4-html/olympic.html Olympic= read.table('olympic.txt', header=TRUE, dec=",") Olympic= Olympic[,1:10] # -------------------------------------------------------- # 2. Pre-traitement # -------------------------------------------------------- # a) p= ncol(Olympic) # nb de variables n= nrow(Olympic) # nb d'individus Mcol= colMeans(Olympic) # Ou en utilisant la commande apply, qui nous permet de calculer la moyenne colonne par colonne : # Mcol= apply(Olympic, 2, mean) # On cree maintenant une matrice a n lignes qui contient Mcol sur chacune de ses lignes : Mcolmat= matrix(Mcol, nrow=n, ncol=p, byrow=TRUE) # On peut alors recentrer les donnees : Oly_centre= Olympic-Mcolmat # b) # Ecart type empirique des colonnes de Oly_centre : norme= apply(Oly_centre, 2, sd) norme # Les variables ne sont pas homogenes. Par exemple, la performance au 100 metre a une variance nettement plus faible que celle au javelot. # ---> On va plutot faire une ACP normee pour comparer des chose comparables norme_mat= matrix(norme, nrow=n, ncol=p, byrow=TRUE) Oly_renorm= Oly_centre/norme_mat # -------------------------------------------------------- # 3. ACP : Representation des athletes # -------------------------------------------------------- # a) acp_olympic= princomp(Oly_renorm, scores=TRUE) acp_olympic # Les standard deviations, ce sont les racines carrees des valeurs propres lambda_j de la matrice de covariance # b) Oly_renorm= as.matrix(Oly_renorm) matcov= (1/n)*(t(Oly_renorm)%*%Oly_renorm) matcov2= (n-1)*cov(Oly_renorm)/n # Ce sont les memes matrices. cov calcule donc (t(Oly_renorm)%*%Oly_renorm)/(n-1) # c) # La commande summary nous donne davantage d'informations : summary(acp_olympic) # En plus des racines des valeurs propres, summary nous donne aussi : # - les variances relatives, cad lambda_j / somme totale des lambda_k (il s'agit de la proportion de variance expliquée par chacune des composantes principales) # - les sommes cumulees des variances relatives, cad sum_{i=1}^q lambda_j / somme totale des lambda_k pour tous les q possibles entre 1 et 10 (il s'agit de la proportion de variance expliquée par le sous-espace vectoriel engendré par les q premières composantes principales) # On va chercher a retrouver les standard deviations calculees par princomp : valpro= eigen(matcov)$values # Calcul des valeurs propres de la matrice de covariance, cad les lambda_j sqrt(valpro) # On retrouve bien les memes valeurs # On peut egalement retrouver les variances relatives avec : (acp_olympic$sdev**2)/sum((acp_olympic$sdev**2)) # Et les somme cumulees des variances relatives avec : cumsum((acp_olympic$sdev**2)/sum((acp_olympic$sdev**2))) # Une visualisation classique realisee apres une ACP est celle dite du "screen plot" : par(mfrow=c(1,2)) plot(1:10, (acp_olympic$sdev**2)/sum((acp_olympic$sdev**2)), xaxt="n", type='b', ylim=c(0,1), main="part de variance expliquee par les composantes", xlab="composante", ylab="part de variance") axis(1,1:10) plot(1:10, cumsum((acp_olympic$sdev**2)/sum((acp_olympic$sdev**2))), xaxt="n", type='b', ylim=c(0,1), main="cumul des parts de variance expliquee par les composantes", xlab="composante", ylab=" cumul des parts de variance") axis(1,1:10) # Permet d'observer l'influence du nombre de composantes sur la variance expliquee (cad la qualite du modele) # Ici, on voit que les trois premieres composantes suffisent a expliquer pres de 80% de la variance # d) acp_olympic$loadings # Les loadings sont les coordonnees des vecteurs propres de la matrice de covariance, en fonction des variables initiales # On peut retrouver ces vecteurs propres avec la commande : vecpro= eigen(matcov)$vectors vecpro # NB : acp_olympics$loadings n'affiche que les coordonnees "franchement non nulles" pour pouvoir interpreter les axes, ici peu parlant car presque tout significatifs, mais si on a plein de rien et juste 0.5 pour le 100 metre et -0.5 pour le disque, ca veut dire que cet axe est en gros la difference entre la performance au 100 metre et la performance du disque # Pour voir tous les non nuls vraiment on tape, par exemple pour la composante 10, acp_olympic$loadings[,10] # e) # L'ajout de l'option scores=TRUE nous permet de recuperer les "scores", cad les coordonnees des individus dans les nouvelles composantes acp_olympic$scores # Pour verifier qu'on obtient bien les bonnes coordonnees, on peut se souvenir # que la premiere colonne est censee etre le produit scalaire pour chaque individu (colonne de taille n) # des p variables de cet individu par le vecteur propre de la composante 1. # Par exemple, pour l'individu 1 sur la 1ere composante principale : # Calcul du produit scalaire sum(Oly_renorm[1,]*acp_olympic$loadings[,1]) # A comparer avec acp_olympic$scores[1,1] # Pour l'individu 2 sur la 1ere composante principale : sum(Oly_renorm[2,]*acp_olympic$loadings[,1]) # A comparer avec acp_olympic$scores[2,1] # On obtient bien le meme resultat a chaque fois. # f) plot(acp_olympic$scores[,1], acp_olympic$scores[,2], type='n', main="Position des individus dans le premier plan factoriel", xlab="1ere composante", ylab="2eme composante") abline(h=0, v=0) text(acp_olympic$scores[,1], acp_olympic$scores[,2], labels=rownames(Olympic)) # NB : le premier axe est engendre principalement (cf la sortie de loadings) par le 110 m haies et le 100 metre en positifs. # Les plus mauvais temps, ie les plus lents a ces deux epreuves sont bien les concurrents 17, 30, 31, 32, 33, 28... # -------------------------------------------------------- # 4. ACP : Representation des epreuves # -------------------------------------------------------- # a) #le vecteur est le vecteur dont la composante i est le coefficient de corrélation linéaire entre les observations de la variable initiale n°1 et les observations sur chacune des nouvelles variables créées. #le vecteur ainsi créé est de norme 1. # b) Cette fois, c'est la corrélation entre les observations sur la première composante principale et chacune des variables initiales. #c) c1<-cor(Oly_renorm,acp_olympic$scores[,1]) c2<-cor(Oly_renorm,acp_olympic$scores[,2]) plot(c1,c2,xlim=c(-1,+1),ylim=c(-1,+1),type="n") abline(h=0,v=0) text(c1,c2,labels=colnames(Olympic),cex=0.5) xt=cos(2*pi*seq(0,1,0.001)) yt=sin(2*pi*seq(0,1,0.001)) points(xt,yt) # -------------------------------------------------------- # 5. Vue d'ensemble # -------------------------------------------------------- # a) biplot(acp_olympic) # En plus de la representation des individus dans le premier plan factoriel, biplot affiche les loadings des anciennes variables # dans les deux premieres composantes. Sur ces deux composantes, on voit des groupes d'epreuves apparaitrent : # - le 100m et le 110m ont des loadings tres similaires (epreuves de sprint) # - le javelot, le disque et le poids sont tres proches (epreuves de lancer) # - la perche, la hauteut et la longueurs sont egalement similaires (epreuves de saut) # L'etude de la matrice de correlation via l'ACP peut donc permettre de faire apparaitre des groupes de variables similaires # b) plot(acp_olympic) #cela trace le diagramme en bâtons des valeurs propres rangées par ordre décroissant. #ce graphique peut permettre d'aider à sélectionner un nombre de composantes à conserver. En effet, en regardant le premier moment où il y a un changement de pente dans la décroissance, cela nous indiquera le nombre de composantes à conserver. # c) #refaire la même chose que précédemment mais à partir de Oly_centre # d) #l'argument cor=TRUE permet de dire au logiciel de procéder à une analyse en composantes principales à partir des données initiales, ce qui permet de ne pas procéder au pré-traitement. Si on ne met pas cor=TRUE, par défaut cela sera une analyse en composantes principales non normées qui sera faite, autrement dit il y aura un pré-traitement de centrage des données, mais pas de normalisation. #Autrement dit, il n'est pas nécessaire de faire les pré-traitements en amont, la fonction princomp les fait toute seule. A vous de spécifier si vous voulez une analyse non normée ou normée (cor=TRUE).