#################### TD ACP exercice 1,  2018 ###
########## UE 7 (Statistiques), Master SVS    ###
########## Patrick Coquillard, 2018           ###
#-----------------------------------------------#


# Partie I : Composante principale et régression linéaire -----

library(MASS)
# mvrnorm(n = 1, mu, Sigma, tol = 1e-6, empirical = FALSE, EISPACK = FALSE)
x <- mvrnorm(1000, c(0,0), matrix(c(5,1.5,1.5,2), 2,2))
x
plot(x, xlim=c(-8,8), ylim=c(-8,8), xlab="X", ylab="Y")

# maintenant les deux régressions linéaires (y = f(x) et x = f(y))
lin.model = lm(x[,1] ~ x[,2])
abline(0, lin.model$coefficients[2], col="blue")
lin.model = lm(x[,2] ~ x[,1])
abline(0, lin.model$coefficients[2], col="blue")

# on fait une ACP pour avoir la composante principale
pca = prcomp(x)
# on récupère la matrice de covariance
r = pca$rotation

# on trace au moyen du coef beta de la régression  dans la matrice de 
# covariance (r) que l'on obtient au moyen de son premier vecteur 
# propre : beta = Vy/Vx

abline(0, r[2,1]/r[1,1], col="red")

# Note -----
# Notez bien que la pente (ß) de la composante n'est pas la moyenne arithmétique 
# des pentes des régressions linéaires, mais la moyenne géométrique (mg)
# des pentes, c.à.d. : mg = sqrt(a*b). 
# On a effectué, en fait, une régression dite "orthogonale" au moyen des 
# "carrés totaux" et non une regression "linéaire" au moyen des 
# "carrés ordinaires" fournis par la fonction lm() habituelle.

# L'intercept de la régression orthogonale s'obtient simplement par (cas
# général) :  ß0 = y¯- ßx¯

# Et si on est dans un espace > 2 ? Et bien c'est la même chose !
# r <- prcomp(cbind(x,y))$rotation
# ß <- -r[-ncol(r),ncol(r)] / r[ncol(r),ncol(r)]
# ß contient toutes les pentes de y = ß1x1 + … + ßkxk.
# A nouveau on a : ß0 = y¯ - x¯ß



# Partie II : une ACP "de base !"-------------------------

#Lecture et preparation des donnees
temperatures <- read.table(file.choose(),sep="\t",header =T)
summary(temperatures)
head(temperatures)     
                                                       
# Maintenant on centre et reduit : x = ((x-mu)/sigma)
t.cr <- scale(temperatures[2:17], center=TRUE, scale=TRUE) 
t.cr
mois <- t.cr[,c(1:12)]      # on oublie les dernieres colonnes
mois
pairs(mois)              # equivalent d'un plot
cov(mois)                # matrice de var-covar
cor(mois)                # matrice des correlations


# On proccède à l'ACP avec le package "stats" standard
acp <- prcomp(mois,scale. = T)# calcule les variances avec (N-1) pour effectifs 
attributes(acp)          # liste du contenu d'acp
summary(acp)

#histogramme des variances = lambdas)
plot(acp, xlab = "Composantes principales (PCi)") 

x <- acp$x ; x           # coordonnées factorielles des objets (après rotation)
plot(x[,1],x[,2], main="position des villes" )#Plot des villes dans le plan 1-2

cov(x)                   # La diagonale des coord fact
                         # contient les variances (lambda)
                         # cov(x) est la matrice diagonale de acp$sdev^2
diag(acp$sdev^2)         # meme resultat...
#par exemple 3.0954² =(SD de PC1)² = 9.5815
                         
round(cov(x))            # l'arrondi montre les axes importants (deux premiers)
                         # et l'indépendance des composantes principales 
 
# tracé des corrélation "à la main" -------------------------------
a <- seq(0,2*pi,length=100)
plot(cos(a),sin(a),type='l',lty = 3, xlab = 'comp 1', ylab = 'comp 2',
      main = "Cercle des corrélations")

r <- acp$r;r                # matrice des vecteurs propres
v <- t(acp$r)               # on transpose
v <- v[1:2,];v              # on récupère les deux premières lignes

# acp$sdev contient les racines carrées des valeurs propres = ecart-types
arrows(0, 0, acp$sdev[1]*v[1,], acp$sdev[2]*v[2,], col = 'red')
text(acp$sdev[1]*v[1,], acp$sdev[2]*v[2,], labels=colnames(v))

# tracé automatique ------------------------------------------------   
require(graphics)
biplot(acp, expand = 1, scale = T)


# Partie III : Autre possibilité ----------------------------------
# verifier que FactoMineR est installee
# sinon  install.packages("FactoMineR", dependencies = TRUE)
#------------------------------------------------------------------#
library(FactoMineR)                                                                                                                           
PCA(temperatures[,2:17])

# a comparer avec  :
acp <- PCA(temperatures[,2:17], quanti.sup= c(13,14,15,16))
acp <- PCA(temperatures[,2:17], quanti.sup= c(13), axes = c(1,3) )
attributes(acp)
acp$call
acp$var$cos2

#------------- Un deuxième jeu de donnees--------------------------

temperature <- read.table("http://factominer.free.fr/livre/temperat.csv",
               header=TRUE, sep=";", dec=".", row.names=1)
res <- PCA(temperature, ind.sup = 24:35, quanti.sup = 13:16, quali.sup = 17)
# !! verifier ce que representent ind.sup, quanti.sup, quali.sup
# Noter ce que contient res : 
res

# habillage défini la colonne qualit. qui défini la coloration d'affichage !
plot.PCA(res, choix="ind", habillage = 17)

dimdesc(res)
res$eig
# resultats sur les individus
res$ind
res$ind.sup

# resultats sur les descripteurs
res$var
res$quanti.sup
res$quali.sup

# centre et reduit les donnees (17 premieres villes)
scale(temperature[1:23,1:16])

#matrice de correlation
cor(temperature[1:23,1:16])

# on trace maintenant les ellipses de chacun des 4 cas
# autour des barycentres et de diamètres égaux aux intervalles
# de confiance des projections sur les deux premiers axes (par défaut).
#seules les villes actives sont prises en compte !

concat.data <- cbind.data.frame(temperature[1:23,17],res$ind$coord)
ellipse.coord <- coord.ellipse(concat.data, bary=TRUE)

# cex est un parametre fixant la taille des caracteres 
plot.PCA(res, habillage =17, ellipse=ellipse.coord, cex=0.8)

#------------- The End ----------------------#