######################################################################
# TP 1 : Analyse en Composantes Principales (ACP) avec le logiciel R 
#       (code R avec commentaires)
######################################################################


#===================================
# Jeu de donnees : eaux petillantes
#===================================
# Chargement du jeu de donnees
eaux <- read.table("../data/eaux.txt", header=TRUE,row.names=1,sep=";")
eaux
##          intensite.emission.de.bulles saveur.salee appreciation.globale
## St Yorre                          3.9          6.4                  2.9
## Vichy                             1.4          6.0                  2.8
## Quezac                            5.1          4.7                  3.5
## Salvetat                          2.9          4.1                  3.4
## Perrier                           8.2          4.9                  2.8
# Chargement du package FactoMineR
library(FactoMineR) # permet de charger le package "PCAmixdata" 
                    # afin de pouvoir l'utiliser par la suite
help(PCA)
# Mise en oeuvre de l'ACP
res <-PCA(eaux)    # tous les calculs de l'ACP sont stockes dans l'objet "res"

                  # NB : par défaut les graphiques des plans factoriels 1-2 
                  # sont affichés a l'écran : 
                  #   - projection des individus sur le plan 1-2
                  #   - cercle des corrélations du plan 1-2
                  
res <- PCA(eaux,graph=FALSE) # idem sans les graphiques 
      
res # permet de voir la liste de l'ensemble des sorties numeriques disponibles
## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 5 individuals, described by 3 variables
## *The results are available in the following objects:
## 
##    name               description                          
## 1  "$eig"             "eigenvalues"                        
## 2  "$var"             "results for the variables"          
## 3  "$var$coord"       "coord. for the variables"           
## 4  "$var$cor"         "correlations variables - dimensions"
## 5  "$var$cos2"        "cos2 for the variables"             
## 6  "$var$contrib"     "contributions of the variables"     
## 7  "$ind"             "results for the individuals"        
## 8  "$ind$coord"       "coord. for the individuals"         
## 9  "$ind$cos2"        "cos2 for the individuals"           
## 10 "$ind$contrib"     "contributions of the individuals"   
## 11 "$call"            "summary statistics"                 
## 12 "$call$centre"     "mean of the variables"              
## 13 "$call$ecart.type" "standard error of the variables"    
## 14 "$call$row.w"      "weights for the individuals"        
## 15 "$call$col.w"      "weights for the variables"
# COMMENTAIRES : par défaut, c'est une ACP centrée réduite qui est faite avec 
# le même poids (1/n) pour tous les individus de l'étude (ici les eux minérales 
# pétillantes). L'espace des individus est R^3 car p=3 variables quantitatives 
# ont été mesurées sur chaque individu. 
# L'espace des variables est R^5 car chaque variable a été mesurée sur n=5 
# individus.
# NB : il s'agit d'un tout petit jeu de données (n=5, p=3) qui contient donc 
# seulement 15 valeurs numériques. Cependant, on va voir que l'ACP permet de 
# faire ressortir très rapidement et visuellement les grandes informations 
# contenues dans ce jeu de données.
# Choix du nombre d'axes à retenir
round(res$eig,digit=2) 
##        eigenvalue percentage of variance cumulative percentage of variance
## comp 1       1.73                  57.78                             57.78
## comp 2       1.06                  35.45                             93.23
## comp 3       0.20                   6.77                            100.00
                       # Permet d'afficher les valeurs propres, 
                       # les pourcentages de variances expliquées par 
                       # axe, puis cumulés.
      
# Graphique de l'ébouli des valeurs propres
barplot(res$eig[,1], main="Eigenvalues", names.arg=1:nrow(res$eig)) 
abline(h=1,col=2,lwd=2)

# COMMENTAIRES : comme on est en ACP normée, on peut utiliser le critère 
# de Kaiser qui dit que seuls les axes factoriels associés à une valeur propre 
# (eigenvalue) plus grande que 1 sont des axes intéressants à retenir.
# En effet, en ACP normée, les variables ont été réduites (divisées  par leur
# écart-type), ainsi leur "nouvelle" variance est 1. 
# Les variables synthétiques (composantes principales) construites par l'ACP 
# ne seront intéressantes que si elles ont une dispersion plus importante
# que celle de n'importe quelle variables de l'étude. Or par construction, 
# la variance d'une composante principale est égale à la valeur propre de l'axe 
# correspondant. D'où le critère de Kaiser qui préconise de ne retenir que 
# les axes ayant une valeur propre supérieure à 1.
#
# NB 1 : si on a une valeur propre inférieure à 1, mais proche de 1, il est 
# recommandé de regarder tout de même quelle est l'information apportée par 
# l'axe correspondant.
#
# NB 2 : si on une valeur propre plus grande que 1 (et proche de 1) qui 
# correspond à un axe factoriel pour lequel on n'arrive pas à faire des 
# commentaires intéressants, il est recommandé de ne pas retenir  cet axe 
# factoriel pour la suite des interprétations.
#
# NB 3 : il y a au maximum p valeurs propres car on peut fabriquer au maximum 
# p composantes principale (=nouvelles variables synthétiques, deux à deux non 
# corrélées, combinaisons linéaires des variables initiales) à partir de p 
# variables initiales.
#
# Ici on obtient les sorties numériques suivantes :
#       Eigenvalue Proportion Cumulative
# dim 1       1.73      57.78      57.78
# dim 2       1.06      35.45      93.23
# dim 3       0.20       6.77     100.00
#
# On a donc les deux premières valeurs propres (les valeurs propres sont 
# toujours rangées par ordre décroissant) supérieures à 1. Selon le critère 
# de Kaiser, il convient donc de retenir 2 axes factoriels. Le premier axe 
# factoriel permet d'expliquer 57,78% de l'inertie (ou de la variance), le 
# second axe factoriel permet d'expliquer 35,45% d'inertie supplémentaire. 
# Ainsi, en considérant le plan 1-2, on récupère 93,23% de l'information 
# (c'est à dire en quoi les individus diffèrent = inertie avec un regard 
# d'algèbre linéaire = variance avec un regard de statisticien).
#
# NB 4 : la troisième valeur propre est très inférieure à 1. Ce troisième axe 
# factoriel ne porte donc pas d'information  pertinente, les 6.77% 
# d'inertie expliquée correspondent plutôt à du bruit (fluctuation aléatoire).
# Graphiques des individus et des variables sur le plan factoriel 1-2
?plot.PCA # permet d'afficher l'aide de la commande "plot.PCA"

plot(res,axes=c(1,2),choix="var") # on retrouve ici le cercle des corrélations (plan 1-2)

# Commentaires du cercle des corrélations du plan 1-2 :
#
# Il s'agit de la projection des variables (qui sont de "taille"longeur" 1 dans l'espace 
# R^5 auquel elles appartiennent) sur le sous-espace de dimension 2 de R^5 qui 
# apporte le plus d'information.
#
# 1) Les 3 variables sont bien représentées dans ce plan 1-2 (ce qui n'est 
# pas une surprise car il contient 93,23% de l'information) car leurs
# projections sont proches de la circonférence du cercle.
#
# 2) Sur l'axe 1 (horizontal), on observe une opposition entre la variable 
# "saveur.salee" (projetée à gauche) et la variable "appreciation.globale" 
# (projetée à droite). Cela signifie que ces deux variables sont fortement 
# corrélées linéairement et négativement entre elles : plus la saveur salée
# est forte, moins l'eau pétillante est appréciée globalement ; a contrario, 
# les eaux les plus appréciées seront celles qui ont une faible saveur salée.
#
# 3) L'axe 2 (vertical) est essentiellement caractérisé par la variable 
# "intensite.emission.de.bulles" qui est projetée vers le haut. 
# Cette variable  (en projection) est orthogonale aux deux autres variables. 
# Cela signifie qu'il n'y a pas a priori des liaisons linéaires entre 
# la variable "intensite.emission.de.bulles" et les deux variables 
# "saveur.salee" et "appreciation.globale". Ainsi, l'intensité d'émission 
# des bulles n'a pas a priori d'impact sur l'appréciation globale de l'eau 
# minérale (ni sur la saveur salée).
#
# Il convient maintenant de projeter les individus (eaux pétillantes) sur le 
# plan factoriel 1-2 associé au sous-espace de dimension 2 de R^3 
# (espace des individus) qui apporte le  plus d'information (maximisation de
# l'inertie du nuage projeté). Ceci nous permettra alors de caractériser 
# les individus en fonction des variables.

plot(res,axes=c(1,2),choix="ind") # on retrouve ici le graphique des individus (plan 1-2)

# Commentaires du plan factoriel 1-2 des individus :
#
# 1) Sur l'axe 1 (horizontal), on voit une opposition très claire entre les 
# aux pétillantes, à gauche, "St Yorre" et "Vichy", et, à droite, "Quezac" et 
# "Salvetat". Au vu cercle des corrélations du plan 1-2, on en déduite que
#      - "St Yorre" et "Vichy" sont des eaux qui ont une saveur salée 
#         importantes et qui n'ont pas été globalement appréciées. 
#      - C'est l'inverse pour "Quezac" et "Salvetat" qui sont des eaux 
#         pétillantes qui ont été globalement bien appréciées et qui n'ont pas 
#         une forte saveur salée.
#      - L'eau pétillante "Perrier" se projette au centre (sur l'axe 
#         horizontal), il s'agit donc d'une eau qui a une saveur salée moyenne 
#         et dont l'appréciation globale a été aussi été moyenne 
#         (moyenne = par rapport aux  cinq eaux pétillantes présentes dans
#          l'étude).
#        
# 2) Sur l'axe 2 (vertical), on va pouvoir caractériser les eaux pétillantes 
# par rapport à leur intensité d'émission  de bulles.
#       - "Perrier" (en haut) diffère très clairement des quatre autres eaux 
#          pétillantes (en bas) : cela signifie que la caractéristique majeure 
#          de "Perrier" est que c'est l'eau pétillante ayant une très forte 
#          intensité d'émission de bulles (par rapport aux 4 autres eaux 
#          pétillantes).
#       - Les eaux "Vichy" et "Salvetat" (en bas) sont celles présentant 
#         la moindre intensité d'émission de bulles.
#       - "St Yorre" et "Quezac" sont dans la moyenne concernant cette 
#         caractéristique d'émission de bulles.


# Au delà des sorties graphiques, des sorties numériques sont disponibles 
# et peuvent être utiles pour quantifier certaines observations graphiques.
# Sorties numeriques pour les individus et pour les variables
# Remarque : par défaut les calculs sont faits sur les 5 premiers axes. 
# (NB : dans cette étude, on peut fabriquer au maximum que 3 axes factoriels).

res$ind # permet d'afficher les sorties numériques associées aux individus : coordonnées, contributions, cosinus carrés.
## $coord
##                 Dim.1       Dim.2      Dim.3
## St Yorre -1.422259675 -0.06130521  0.5185196
## Vichy    -1.520947213 -0.90545543 -0.3700892
## Quezac    1.422493513 -0.13598570  0.5801695
## Salvetat  1.518799529 -0.83922841 -0.4280543
## Perrier   0.001913846  1.94197475 -0.3005456
## 
## $cos2
##                 Dim.1       Dim.2      Dim.3
## St Yorre 8.812339e-01 0.001637300 0.11712881
## Vichy    7.074044e-01 0.250711156 0.04188440
## Quezac   8.507138e-01 0.007774444 0.14151179
## Salvetat 7.221493e-01 0.220488754 0.05736194
## Perrier  9.485212e-07 0.976607789 0.02339126
## 
## $contrib
##                 Dim.1       Dim.2     Dim.3
## St Yorre 2.334112e+01  0.07067623 26.463285
## Vichy    2.669268e+01 15.41745729 13.481126
## Quezac   2.334880e+01  0.34774834 33.130132
## Salvetat 2.661735e+01 13.24460325 18.034795
## Perrier  4.226472e-05 70.91951488  8.890661
## 
## $dist
## St Yorre    Vichy   Quezac Salvetat  Perrier 
## 1.515072 1.808341 1.542263 1.787257 1.965095
# NB : les coordonnées peuvent être utiles dans le cas où deux individus sont projetés 
# au même endroit sur le plan factoriel et qu'il est difficile de lire leurs noms graphiquement.

round(res$ind$cos2,digits=3) # uniquement les cosinus carres (arrondis à 3 chiffres après la virgule)
##          Dim.1 Dim.2 Dim.3
## St Yorre 0.881 0.002 0.117
## Vichy    0.707 0.251 0.042
## Quezac   0.851 0.008 0.142
## Salvetat 0.722 0.220 0.057
## Perrier  0.000 0.977 0.023
# Commentaires : les cosinus carrés nous donnent la qualité de représentation (en %) 
# des individus  sur les axes factoriels.
# Par exemple : 70.7 % de l'information de "Vichy" est recupérée sur l'axe 1, 25,1% sur l'axe 2.
# Vu que les axes sont orhtogonaux, on a 70,7+25,1=95,8% de l'information de "Vichy" 
# qui est fournie par le plan 1-2.
# "Perrier" est extrêmement bien représentée sur l'axe 2 (97,7%).
# On observe ici que les n=5 individus sont très bien représentés sur le plan 1-2.
#  Ceci n'est pas surprenant car la qualité globale de plan 1-2 est de 93,23% 
# (d'inertie ou variance expliquée).
#
# NB : Ne jamais faire des commentaires sur des individus qui seraient mal représentés 
# sur le plan factoriel que l'on regarde car on ne ferait que des fausses déductions... 
 

round(res$ind$contrib, digits=1) # uniquement les contributions en % (arrondies à un chiffre après la virgule)
##          Dim.1 Dim.2 Dim.3
## St Yorre  23.3   0.1  26.5
## Vichy     26.7  15.4  13.5
## Quezac    23.3   0.3  33.1
## Salvetat  26.6  13.2  18.0
## Perrier    0.0  70.9   8.9
# Commentaires : On voit que 4 individus contribuent formement à l'axe 1 
# "St Yorre", "Vichy", "Quezac" et "Salvetat" (avec des contributions de 
# l'ordre de 25%). "Pierrier" n'a pas du tout contribué à la fabrication de 
# ce premier axe factoriel.
# Pour l'axe 2, c'est essentiellement "Perrier" qui a fortement contribué 
# à sa fabrication (70,9%) et dans une moindre mesure "Vichy" et "Salvetat" 
# (avec des contributions de l'ordre de 15%). 
#
# NB : je rappelle que l'on ne souhaite pas que seul un individu (ou peu 
# d'invidus lorsque n est grand) contribue à la fabrication du (ou des) 
# premier(s) axe(s) factoriels car cela rendrait alors les résultats de l'étude 
# instables dans le sens où la présence ou pas de cet individu dans  
# les données va entrainer une modification dans la construction des premiers axes factoriels.
# Ici, ce n'est pas le cas. On a 4 individus sur 5 (80% des individus) qui 
# contribuent au premier axe factoriel, et 3 sur 5 (60%) qui contribuent au
# second axe. C'est normal que "Perrier" qui n'a pas contribué à la fabrication 
# de l'axe 1 (qui contient 57,78% de l'information) ait tendance à contribuer 
# fortement à la constitution de l'axe 2 (qui contient 35,45% de l'information).


res$var # permet d'afficher les sorties numériques associées aux variables :  coordonnées, contributions, cosinus carrés.
## $coord
##                                   Dim.1       Dim.2     Dim.3
## intensite.emission.de.bulles  0.2644680  0.95324075 0.1462489
## saveur.salee                 -0.9462657 -0.08864987 0.3110022
## appreciation.globale          0.8763031 -0.38341534 0.2916943
## 
## $cor
##                                   Dim.1       Dim.2     Dim.3
## intensite.emission.de.bulles  0.2644680  0.95324075 0.1462489
## saveur.salee                 -0.9462657 -0.08864987 0.3110022
## appreciation.globale          0.8763031 -0.38341534 0.2916943
## 
## $cos2
##                                   Dim.1       Dim.2      Dim.3
## intensite.emission.de.bulles 0.06994335 0.908667926 0.02138873
## saveur.salee                 0.89541882 0.007858799 0.09672238
## appreciation.globale         0.76790713 0.147007322 0.08508555
## 
## $contrib
##                                  Dim.1      Dim.2    Dim.3
## intensite.emission.de.bulles  4.035342 85.4385366 10.52612
## saveur.salee                 51.660687  0.7389326 47.60038
## appreciation.globale         44.303971 13.8225308 41.87350
# NB : les coordonnées peuvent être utiles dans le cas où deux 
# variables sont projetées au même endroit sur le cercle des 
# corrélations et qu'il est difficile de lire leurs noms graphiquement.
             
round(res$var$cos2,digits=3) # uniquement les cosinus carres (arrondis à 3 chiffres après la virgule)
##                              Dim.1 Dim.2 Dim.3
## intensite.emission.de.bulles 0.070 0.909 0.021
## saveur.salee                 0.895 0.008 0.097
## appreciation.globale         0.768 0.147 0.085
# Comme pour les individus, les cosinus carrés nous donnent la qualité (en %) 
# de réprésentation des variables sur les axes factoriels.
# Par exemple, la variable "intensite.emission.de.bulles" est très mal 
# représentée sur l'axe 1 (7% seulement de l'information de cette variable est  
# portée par cet axe), mais elle est très bien représentée sur l'axe 2 (90,9% 
# de l'information de cette variable est fournie par l'axe vertical).
# Globalement sur le plan principal (1-2), cette variable 
# "intensite.emission.de.bulles" a une qualité de représentation égale à 
# 7% + 90,9% = 97,9%, elle est donc quasiment parfaitement représentée sur 
# le plan 1-2.
#
# NB : on peut sommer les qualités de représentation de l'axe 1 et de l'axe 2 
# pour avoir la qualité de représentations ur le plan 1-2, car les axes sont
# construits orthogonaux entre eux, et il n'y a donc pas de redondance 
# d'information entre les axes factoriels.
#
# Les deux autres variables "saveur.salee" et "appreciation.globale" sont 
# très bien représentée sur l'axe 1 (resp. 89.5% et 76,8%) et mal représentée 
# sur l'axe 2 (resp. 0.8% et 14,7%). 

round(res$var$contrib,digits=1) # uniquement les contributions en % (arrondies à un chiffre après la virgule)
##                              Dim.1 Dim.2 Dim.3
## intensite.emission.de.bulles   4.0  85.4  10.5
## saveur.salee                  51.7   0.7  47.6
## appreciation.globale          44.3  13.8  41.9
# Les contributions (en %) des variables à la fabrication des axes factoriels 
# nous indiquent naturellement que seules les variables "saveur.salee" et 
# "appreciation.globale" ont beaucoup contribué à la fabrication de 
# l'axe 1 (resp. 51.6% et 44,3%). La variable "intensite.emission.de.bulles" 
# est la plus forte contributrice à la fabrication de l'axe 2
# (avec un pourcentage de contribution de l'ordre de 85,4%).