# petites explorations autour de l'ACP avec R # Time-stamp: <2007-07-30 20:31:09 chl> # # Note: # voir article Rnews Décembre 2003, "Dimensional Reduction # for Data Mapping", pp. 2-7 (Rnews_2003-3.pdf) # # génération des données utilisées dans l'article: helix <- function(size=30) { # input : size, sample size # output: 3d helix t <- 0:(size-1) z <- sqrt(2)/2*t y <- sin(z) x <- cos(z) cbind(x,y,z) } gauss54 <- function(s=15, sig=0.02) { # input : s, sample size (per Gaussian) # sigma, cluster_covariance_ # output: 5 4d Gaussian clusters require(MASS) simp4d <- matrix(c(0,0,0,0, 1,0,0,0, 1/2,sqrt(3)/2,0,0, 1/2,sqrt(3)/6,sqrt(2/3),0, 1/2,sqrt(3)/6,sqrt(2)/(4*sqrt(3)),sqrt(5/8)), 5,4,byrow=T) # simplex4d can easily be checked using 'dist' rbind(mvrnorm(simp4d[1,],S=diag(4)*sig,n=s), mvrnorm(simp4d[2,],S=diag(4)*sig,n=s), mvrnorm(simp4d[3,],S=diag(4)*sig,n=s), mvrnorm(simp4d[4,],S=diag(4)*sig,n=s), mvrnorm(simp4d[5,],S=diag(4)*sig,n=s)) } l9 <- matrix(rep(1:9,9),9,9) h30 <- helix() # g54, last column is class label g54 <- cbind(gauss54(),rep(1:5,e=15)) # iris data is a standard data set data(iris) # check for duplicates else distance algorithms 'blow up' iris.data <- unique(iris[,1:4]) iris.class <- iris[row.names(iris.data),5] # some colours for the maps classcol <- c("red","blue","green","black","yellow") # usage par(mfrow=c(3,3)) out.pca <- prcomp(g54[,1:4]) plot(out.pca$x[,1:2],main="prcomp") dist.data <- dist(g54[,1:4]) out.pcoa <- cmdscale(dist.data) randstart <- matrix(runif(nrow(g54[,1:4])*2),nrow(g54[,1:4]),2) out.sam <- sammon(dist.data,y=randstart,tol=1e-6,niter=500) out.iso <- isoMDS(dist.data,y=randstart,tol=1e-6,maxit=500) plot(out.pcoa,main="cmdscale") plot(out.sam$points,ylab="",xlab="", main=paste("sammon s=",as.character(round(out.sam$stress,4)))) plot(out.iso$points,ylab="",xlab="", main=paste("IsoMDS s=",as.character(round(out.iso$stress,4)))) out.ica <- fastICA(g54[,1:4],2,alg.typ="deflation",fun="logcosh", alpha=1,row.norm=FALSE,maxit=200,tol=0.0001) plot(out.ica$[,1:2]) par(mfrow=c(2,1)) plot(as.matrix(dist.data),as.matrix(dist(out.sam$points)), main="Distorsion plot: Sammon map",ylab="2d distance", xlab="nd distance") plot(as.matrix(dist.data),as.matrix(dist(out.iso$points)), main="Distorsion plot: IsoMDS map",ylab="2d distance", xlab="nd distance") library(class) gr <- somgrid(topo="hexagonal",xdim=5,ydim=5) out.som <- SOM(g54[,1:4],gr, alpha=list(seq(0.05,0,len=1e4),seq(0.02,0,len=1e5)), radii=list(seq(8,1,len=1e4),seq(4,1,len=1e5))) plot(out.som$grid,type="n") symbols(out.som$grid$pts[,1],out.som$grid$pts[,2], circles=rep(0.4,25),inches=FALSE,add=TRUE) bins <- as.numeric(knn1(out.som$codes,g54[,1:4],0:24)) points(out.som$grid$pts[bins,]+rnorm(nrow(g54),0,0.1), col=classcol[rep(1:5,e=15)],pch=20) # la même chose avec les données `iris` ###################################################################### # Un exemple simple : l'ACP pour des données bivariées ###################################################################### # # Exemple fictif 9 ind. x 4 var. # # Source des données: # Baccini & Besse, Data mining -- 1. Exploration statistique, 2004 (p. 24) # http://www.lsp.ups-tlse.fr # a <- read.table("notes.txt",header=T) # résumer l'information numérique summary(a[-1]) # si besoin, on peut centrer les variables a.sc <- scale(a[-1]) # matrice de covariances et de corrélation a.cor <- cor(a[-1]) a.cov <- cov(a[-1]) # représenter graphiquement les liaisons entre variables pairs(a[-1]) # Note: # voir aussi la fonction cpairs {gclus} qui permet de réarranger # les panels en mettant les corrélations les plus élevées près de # la diagonale, ou d'indiquer les corrélations par des couleurs # matrice de dispersion plus "sexy" library(ggplot2) plotmatrix(a[-1]) # matrice de dispersion plus "sexy" (bis) library(lattice) splom(a[-1]) # matrice de dispersion avec distributions univariées library(car) scatterplot.matrix(a[-1]) # matrice de dispersion + valeurs de corrélation library(pgirmess) pairsrp(a[-1]) # même chose avec plus d'informations library(languageR) pairscor.fnc(a[-1]) detach("package:lattice") detach("package:ggplot2") detach("package:pgirmess") detach("package:languageR") ###################################################################### # ACP à la main (sans utiliser les packages R) # # le but est plutôt pédagogique : on décompose l'analyse étape par # étape en commençant par l'analyse portant sur les variables, puis # les caractéristiques des individus dans le nouveau repère # acp.summary <- function(m, type="cov") { if (type == "cov") x <- cov(m) else if (type == "cor") x <- cor(m) else stop("Unknown option type.") vap <- eigen(x)$values vap.tot <- sum(vap) vap.perc <- vap/vap.tot vep <- eigen(x)$vectors # valeurs propres et vecteurs propres vap.m <- matrix(round(c(vap,vap.perc,cumsum(vap.perc)),3),nrow=length(vap), dimnames=list(paste("F",1:length(vap),sep=""), c("val. pr.","% var.","% var. cum."))) # corrélations variables/facteurs fac.coor <- vep*m varfac <- cor(fac.coor) # nuage des individus ## poids <- 1/dim(m)[1] ## cont.var <- (poids*fac.coor)/vap*100 return(list("Val. prop."=vap.m, "Var. tot."=round(vap.tot,2), "Corr. var./fac."=varfac)) } acp.summary(a[-1]) # Note: # une décomposition en valeurs signulières avec la fonction `svd` # donnera (au signe près) les mêmes vecteurs propres # Illustration dans le cas 2D # Tiré du site de V. Zoonekynd # http://zoonek2.free.fr/UNIX/48_R/05.html N <- 1000 X <- matrix(runif(2*N, -1, 1), nc=2) plot(X) M <- matrix(rnorm(4), nc=2) Y <- X %*% M plot(Y) plot(Y) p <- prcomp(Y)$rotation abline(0, p[2,1] / p[1,1], col="red", lwd=3) abline(0, -p[1,1] / p[2,1], col="red", lwd=3) ###################################################################### # ACP avec princomp ###################################################################### # ACP avec la librairie ADE4