########################################################### # Liste des fonctions sous R 2.14.1 (version ADlin 0.968): # modifié 04/06/2021 labels in PCA plots display for example "Dim 1 (61.26%)" # modifié 16/06/2015 Arrows for the variables in PCA correlations circle # modifié 06/06/2015 MVcut avec l'option graph pour choisir à la main les breaks # modifié 20/05/2015 version en anglais de la fonction acpxqd() pour l'ACP usuelle: la fonction PCA() # modifié 20/05/2015 fonction Hotelling pour détecter T² outliers # modifié 11/12/2012 les return() multiples arguments sont listés # modifié 01/10/2004 fonction mygraphics pour sauvegarder des images en fichiers postscrip, pdf ou jpeg # modifié le 28/09/2004 : acpxqd et afc ont les plots avec background couleur, bgpar=... # modifié le 02/04/2004 : acpxqd et afc (cas AFCM) : carte couleur des individus (lignes) # selon les modalités d'une variable qualitative # modifié le 25/11/2003 : afc avec elements supplementaires # .First STATIS RVop STATIS.plot Dcenter invgene Dproj Dcp Dvar # WDop VQop ACPVI acpxqd afc posmultclas star.graph table.val xy.val # codisj Codisjc discri MVcut myps Profils PCA ########################################################### #.First<-function() # STARTUP FONCTION A N'UTILISER QUE SOUS UNIX #{ #X11("-xrm 'splus*colors: white black red 2 green 2 blue'") #} ########################################################### STATIS<-function(X,D=1) { if(!is.list(X))return() k<-length(X) Y<-list(NULL) WD<-list(NULL) n<-nrow(X[[1]]) for(i in 1:k) { Y[[i]]<-Dcenter(X[[i]],D=D)$X WD[[i]]<-WDop(Y[[i]],D=D) } R<-matrix(1,k,k) for(i in 1:(k-1)) for(j in (i+1):k){ R[i,j]<-RVop(WD[[i]],WD[[j]],D=D) R[j,i]<-R[i,j] } diag<-eigen(R,symmetric=T) valp<-diag$values vect<-diag$vectors alpha<-as.vector(abs(diag$vectors[,1])) WDC<-matrix(0,n,n) for(i in 1:k)WDC<-WDC+alpha[i]*WD[[i]] return(list(R=R,valp=valp,alpha=alpha,vect=vect,WDC=WDC,D=D)) } ################################################ RVop<-function(X,Y,D=1) {if(nrow(X)!=ncol(X))return() if(nrow(Y)!=ncol(Y))return() if(nrow(X)!=nrow(Y))return() X<-as.matrix(X) Y<-as.matrix(Y) if(length(D)==1) DD<-rep(1/nrow(Y),nrow(Y)) else DD<-as.vector(D) A<-as.matrix(DD*X) B<-as.matrix(DD*Y) if(sum(A-t(A))!=0)return() if(sum(B-t(B))!=0)return() num<-sum(X*Y) den<-sum(X*X)*sum(Y*Y) RV<-num/sqrt(den) return(RV) } #################################################################### STATIS.plot<-function(resul) {k<-length(resul$valp) par(ask=T) # hit return to see next plot x<-as.vector(resul$vect[,1]*sqrt(resul$valp[1])) y<-as.vector(resul$vect[,2]*sqrt(resul$valp[2])) plot(c(x,0),c(y,0),xlab="axe 1",ylab="axe 2",type="n") text(x,y) arrows(rep(0,k),rep(0,k),x*0.95,y*0.95,angle=15,length=0.15) diagWDC<-eigen(round(resul$WDC,6),symmetric=T) xx<-diagWDC$vectors[,1]*sqrt(diagWDC$values[1]) yy<-diagWDC$vectors[,2]*sqrt(diagWDC$values[2]) plot(xx,yy,xlab="axe 1",ylab="axe 2",type="n") abline(h=0) abline(v=0) text(xx,yy,dimnames(resul$WDC)[[1]]) } ################################################# invgene<-function(A,eps=1e-06) # inverse generalisee de A { A<-as.matrix(A) # decomposition valeurs singulieres de A et calcul de son inverse dans RR valsin<-svd(A) diago<-valsin$d[valsin$d>eps] # browser() if(length(diago) == 0){RR<-matrix(0,ncol(A),nrow(A)) return(RR)} if(length(diago)==1) RR <- as.matrix(valsin$v[, 1:length(diago)])%*%t(as.matrix(valsin$u[, 1:length(diago)]))/diago else RR<-valsin$v[,1:length(diago)]%*%diag(1/diago)%*%t(valsin$u[,1:length(diago)]) #if(length(RR)==1)RR<-as.numeric(RR) return(RR) } ############################################################## Dproj<-function(X,Y,D=1,eps=1e-06,inv=T) # D-proj de Y sur X # sortie du residu Yres=Y-P Y de Yhat=P Y et des coef. beta # X X { X<-as.matrix(X) Y<-as.matrix(Y) if(inv)R<-solve(Dcp(X,D=D)) else R<-invgene(Dcp(X,D=D),eps=eps) Yhat<-WDop(X,Q=R,D=D)%*%Y Yres<-Y-Yhat beta<-R%*%Dcp(X,Y,D=D) return(list(Yhat=Yhat,Yres=Yres,beta=beta)) } ################################################################ Dcp<-function(X,Y=X,D=1) { # signification de Dcp : Diagonalcrossproduct, i.e., X'DY # D est un vecteur des poids, par defaut 1/nrow(X) Y<-as.matrix(Y) X<-as.matrix(X) if(length(D)==1) DD<-rep(1/nrow(X),nrow(X)) else DD<-as.vector(D) for(i in 1:nrow(X)) X[i,]<-DD[i]*X[i,] V<-crossprod(X,Y) dimnames(V)<-NULL return(V) } ################################################################ Dvar<-function(X,D=1,cor=F) # VALUE: # V: matrice des D-variances de X si cor=F, sinon matrice des D-correlations # U: matrice D-centree de X si cor=F, sinon matrice D-centree reduite # mean: vecteur des D-moyennes # var: vecteur des D-variances { mat<-FALSE if(is.matrix(X)) mat<-TRUE X<-as.matrix(X) centrage<-Dcenter(X,D=D) Y<-centrage$X mean<-centrage$mx V<-VQop(Y,D=D) var<-diag(V) if(cor==F){ U<-Y } else { ect<-sqrt(var) U<-sweep(Y,2,ect,FUN="/") V<-VQop(U,D=D) } if(mat){ dimnames(V)<-list(dimnames(X)[[2]],dimnames(X)[[2]]) dimnames(U)<-dimnames(X) } return(list(V=V,U=U,mean=mean,var=var)) } ################################################################# WDop<-function(X,Y=X,Q=1,D=1) { # X,Y vecteurs ou matrices n x p, par defaut Y=X; # operateur des produits scalaires entre individus de X et Y, XQY'D; # Q est la metrique de l'espace des individus (par defaut Q=Ipp);Q matrice # D est la metrique de l'espace des variables de Y:(par defaut 1/n Inn), # D vecteur X<-as.matrix(X) Y<-as.matrix(Y) if(length(D)==1) DD<-rep(1/nrow(Y),nrow(Y)) else DD<-as.vector(D) if(!is.matrix(Q)) W<-X%*%t(Y) else W<-X%*%Q%*%t(Y) for(i in 1:ncol(W)) W[,i]<-DD[i]*W[,i] return(W) } ################################################################## VQop<-function(X,Y=X,Q=1,D=1) { # operateur crossproduct generalise, i.e., X'DYQ; # Q matrice ou vecteur est une metrique sur l'espace des variables de Y, # par defaut Q=1 donne l'identite # si Q=k alors Q=kI # si Q=vecteur alors Q=diag(vecteur); # D est un vecteur, par defaut D=1 donne le vecteur 1/nrow(X); Y<-as.matrix(Y) X<-as.matrix(X) if(length(D)==1) DD<-rep(1/nrow(X),nrow(X)) else DD<-as.vector(D) for(i in 1:nrow(X)) X[i,]<-DD[i]*X[i,] if(! is.matrix(Q)) { QQ<-as.vector(Q) if(length(QQ)==1 & QQ[1]==1){ V<-crossprod(X,Y) return(V) } if(length(QQ)==1 & QQ[1]!=1) QQ<-rep(QQ,ncol(Y)) for(i in 1:ncol(Y)) Y[,i]<-Y[,i]*QQ[i] V<-crossprod(X,Y) return(V) } V<-crossprod(X,Y)%*%Q return(V) } ################################################################### Dcenter<-function(X,D=1) #D-centrage de la matrice X { X<-as.matrix(X) if(all(!is.na(X))){ if(length(D)==1) mx<-apply(X,2,"mean") else {idx<-D*X mx<-apply(idx,2,"sum") } X<-sweep(X,2,mx) } else { mx<-rep(0,ncol(X)) if(length(D)==1){for(i in 1:ncol(X)) mx[i]<-mean(X[,i][!is.na(X[,i])])} else {idx<-D*X mx<-apply(idx,2,"sum") } X<-sweep(X,2,mx) } return(list(X=X,mx=mx)) } ################################################################### Dcentred<-function(X,D=1) { # D-centrage de la matrice X # # Entree # X matrice des variables # D metrique des poids # Sorties # Xc matrice deduite de X, D-centree # Xcr matrice deduite de X, D-centree et reduite # moy vecteur des D-moyennes # var vecteur des D-variances X<-as.matrix(X) n<-nrow(X) p<-ncol(X) n2<-n*n n1<-matrix(1,nrow=n,ncol=1) if(length(D)==1)D<-diag(rep(1/n,n),nrow=n) if(length(D)==n)D<-diag(D,nrow=n) if(length(D)==n2)D<-as.matrix(D) moy<-t(n1)%*%D%*%X Xc<-sweep(X,2,moy) var<-diag(t(Xc)%*%D%*%Xc) ect<-sqrt(var) Xcr<-sweep(Xc,2,ect,FUN="/") return(list(Xc=Xc,Xcr=Xcr,moy=moy,var=var)) } ############################################################## ACPVI<-function(Y, QY = 1, D = 1, X, centrer = T, cor = T, k = 3, impres = T, graph = T) { # # ACPVI de X par rapport a un triplet (Y,QY,D) # # Entrees # # Y matrice des variables a expliquer # QY metrique des u.s. pour Y : si QY=1 alors QY=Idn, si QY est un vecteur alors QY=diag(vecteur), sinon QY=matrice # D metrique des variables, si D=1 alors 1/nIdn, si D est un vecteur alors D=diag(vecteur), sinon D=matrice # X matrice des variables explicatives ou variables instrumemtales # centrer=F on ne centre pas X et Y, sinon centrer=T # cor=F on ne reduit pas, centrage et reduction si cor=T # k nombre de composantes retenues, par defaut 3 # impres impression des resultats si T, si F rien # graph trace des graphiques si T, si F pas de trace # # Sorties # # nomfichX nom du fichier X # nomfichY nom du fichier Y # Xini tableau initial # Yini tableau initial # X tableau apres eventuel centrage et/ou reduction # Y tableau apres eventuel centrage et/ou reduction # Yhat projection de Y sur Im(X) i.e. regression multivariee de Y par X # Yres Y-Yhat residu # QY metrique utilisee pour le triplet de Y # D metrique utilisee pour les u.s. de X et Y # inertiaY inertie du triplet (Y,QY,D) # inertiaACPVI inertie du triplet (X,QX,D) # QX metrique optimale pour X # k rang de la DSD retenue # valp vecteur de toutes les valeurs propres de l'ACPVI # A matrice des k premiers axes de l'ACPVI, des correlations variables de Y composantes si cor=T # C matrice des k premieres composantes principales de l'ACPVI # ZA matrice des correlations de X avec les composantes principales de l'ACPVI si cor=T # # Fonctions exterieures appelees # # Dcentred, invgene, acpxqd # # R.Sabatier # MAJ 03/04/97 # # lecture des donnees nomfichX <- deparse(substitute(X)) if(!is.matrix(X)) stop("X n'est pas une matrice") X <- as.matrix(X) if(length(c(which.inf(X), which.na(X)))) stop("valeurs manquantes ou infinies dans X") Xini <- as.matrix(X) nomfichY <- deparse(substitute(Y)) if(!is.matrix(Y)) stop("Y n'est pas une matrice") Y <- as.matrix(Y) if(length(c(which.inf(Y), which.na(Y)))) stop("valeurs manquantes ou infinies dans Y") Yini <- as.matrix(Y) n <- nrow(X) p <- ncol(X) q <- ncol(Y) if(n != nrow(Y)) return() if(p <= 2) return() if(impres == F) graph <- F if(centrer == F) cor <- F if(cor == T) centrer <- T q2 <- q * q n2 <- n * n # calcul de la metrique QY imet <- as.numeric(1) if(length(QY) == q2) { QY <- as.matrix(QY) imet <- 2 } if(length(QY) == 1) QY <- diag(q) if(length(QY) == q) QY <- diag(QY, nrow = q) # calcul de la metrique D if(length(D) == 1) D <- diag(rep(1/n, n), nrow = n) if(length(D) == n) D <- diag(D, nrow = n) if(length(D) == n2) D <- as.matrix(D) # centrage et/ou reduction centrageX <- Dcentred(X, D = D) meanX <- matrix(centrageX$moy) varX <- matrix(centrageX$var) dimnames(meanX) <- list(dimnames(X)[[2]], "moy") dimnames(varX) <- list(dimnames(X)[[2]], "var") centrageY <- Dcentred(Y, D = D) meanY <- matrix(centrageY$moy) varY <- matrix(centrageY$var) dimnames(meanY) <- list(dimnames(Y)[[2]], "moy") dimnames(varY) <- list(dimnames(Y)[[2]], "var") if(centrer == T) { X <- centrageX$Xc Y <- centrageY$Xc } if(cor == T) { X <- centrageX$Xcr Y <- centrageY$Xcr } # calcul de la metrique optimale QX invVX <- invgene(t(X) %*% D %*% X) VXY <- t(X) %*% D %*% Y QX <- invVX %*% VXY %*% QY %*% t(VXY) %*% invVX # calcul de quelques valeurs optimales associees a l'ACPVI VQY <- t(Y) %*% D %*% Y %*% QY inertiaY <- sum(diag(VQY)) normY <- sum(diag(VQY %*% VQY)) if(impres) { cat("__________________________________________________________________\n" ) cat(" - A.C.P.V.I. de X par rapport a un triplet (Y,QY,D) -\n" ) if(cor == T) { cat(" Pour X : ", nomfichX, "\n") cat(" donnees centrees et reduites\n") cat(" moyenne des variables de X\n") print(t(meanX)) cat(" variance des variables de X\n") print(t(varX)) cat(" Pour Y : ", nomfichY, "\n") cat(" donnees centrees et reduites\n") cat(" moyenne des variables de Y\n") print(t(meanY)) cat(" variance des variables de X\n") print(t(varY)) cat(" inertie du triplet de Y :", format( inertiaY), "\n") } if(centrer == T && cor == F) { cat(" Pour X : ", nomfichX, "\n") cat(" donnees centrees\n") cat(" moyenne des variables de X\n") print(t(meanX)) cat(" variance des variables de X\n") print(t(varX)) cat(" Pour Y : ", nomfichY, "\n") cat(" donnees centrees\n") cat(" moyenne des variables de Y\n") print(t(meanY)) cat(" variance des variables de Y\n") print(t(varY)) cat(" inertie du triplet de Y :", format( inertiaY), "\n") cat("__________________________________________________________________\n" ) } if(centrer == F) { cat(" Pour X et Y donnees non transformees\n") cat(" inertie du triplet de Y :", format( inertiaY), "\n") cat("__________________________________________________________________\n" ) } } # calcul des differentes matrices Yhat <- X %*% invVX %*% VXY # matrice Y apres projection sur Im(X) Yres <- Y - Yhat # matrice residuelle de Y apres projection sur Im(X) varYhat <- matrix(Dcentred(Yhat, D = D)$var) dimnames(varYhat) <- list(dimnames(Y)[[2]], "var.expliquee") if(impres) { cat("__________________________________________________________________\n" ) cat("Variance des variables de Y apres ACPVI\n") print(t(varYhat)) } # ACPVI resACPVI <- acpxqd(Yhat, Q = QY, D = D, centrer = F,k=k, impres = F, graph = F) inertiaACPVI <- resACPVI$inertiaX valp <- resACPVI$valp normYhat <- sum(valp^2) C <- resACPVI$C k <- resACPVI$k A <- resACPVI$A # impression de l'histogramme des valeurs propres if(impres) { cat("__________________________________________________________________\n" ) cat(paste("Inertie de l'ACPVI =", format(inertiaACPVI), "(", format(inertiaACPVI/inertiaY * 100), "%)\n")) cat(paste("Carre de la norme de l'ACPVI =", format(normYhat), "\n")) cat(paste("RV de l'ACPVI =", format(sqrt(normYhat/normY)), "\n")) cat(paste("Carre de la distance entre les deux operateurs =", format(normY - normYhat), "\n")) cat("__________________________________________________________________\n" ) cat("histogramme des valeurs propres (o ou n) ?\n") plth <- scan(quiet=T,quiet=T,"", character(), 1) if(plth == "o") { par(mfrow = c(1, 1)) barplot(valp, space = 2, names = paste("v. p.", 1:q)) title("valeurs propres") } } # choix du nombre d'axes si non donne en argument cat("__________________________________________________________________\n" ) valtab <- matrix(0, q, 4) dimnames(valtab) <- list(format(1:q), c("val.pro.", "% inert.", "% cumul.", " RV")) valtab[, 1] <- round(valp, digits = 4) valtab[, 2] <- round(valp/inertiaACPVI * 100, digits = 2) for(i in 1:q) { valtab[i, 3] <- sum(valtab[1:i, 2]) valtab[i, 4] <- sqrt(sum(valp[1:i]^2)/normY) } if(k == 0) { print(valtab) repeat { cat("Combien d'axes voulez-vous ? (<=", q, ")\n") k <- scan(quiet=T,quiet=T,"", n = 1) if(k != 0) break } } ZA <- t(C) %*% D %*% X dimnames(ZA) <- list(paste("c", format(1:k), sep = ""), dimnames(X)[[2 ]]) dimnames(A) <- list(dimnames(Y)[[2]], paste("a", format(1:k), sep = "") ) dimnames(C) <- list(dimnames(Y)[[1]], paste("c", format(1:k), sep = "") ) # aides a l'interpretations eventuels des variables de Y if(impres == T) { if(imet == 1) { repeat { cat("Aide a l'interpretation des variables (o/n) ?\n" ) aiv <- scan(quiet=T,quiet=T,"", character(), 1) if((length(aiv) == 0) | (aiv == "n")) break else { aivtab0 <- matrix(0, q, 2) aivtab1 <- matrix(0, q, k) aivtab2 <- matrix(0, q, k) dimnames(aivtab0) <- list(dimnames(Y)[[2]], c( "var.exp.", "part")) dimnames(aivtab1) <- list(dimnames(Y)[[2]], paste("CTR", format(1:k), sep = "")) dimnames(aivtab2) <- list(dimnames(Y)[[2]], paste("PART", format(1:k), sep = "")) aivtab0[, 1] <- round(varYhat * 1000, digits = 0) for(j in 1:k) aivtab1[, j] <- round((QY[j, j] * A[, j]^2)/ valp[j] * 1000, digits = 0) for(i in 1:q) aivtab2[i, ] <- round(A[i, ]^2/varYhat[i] * 1000, digits = 0) for(i in 1:q) aivtab0[i, 2] <- sum(aivtab2[i, ]) aivtab <- data.frame(aivtab0, aivtab1, aivtab2) print(aivtab) } } } } # trace des eventuels graphique if(graph == T) { repeat { cat("graphique pour les composantes (o/n) ?\n") pltc <- scan(quiet=T,"", character(), 1) if((length(pltc) == 0) | (pltc == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltch <- scan(quiet=T,"", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltcv <- scan(quiet=T,"", numeric(), 1) par(mfrow = c(1, 1), pty = "s") axespar <- c(pltch, pltcv) plot(C[, axespar], xlab = paste("c", axespar[1], " ", round(valp[axespar[1]], digits = 4), "(", round(valp[axespar[1]]/inertiaACPVI * 100, digits = 2), "%)"), ylab = paste("c", axespar[ 2], " ", round(valp[axespar[2]], digits = 4), "(", round(valp[axespar[2]]/inertiaACPVI * 100, digits = 2), "%)"), type = "n") abline(h = 0) abline(v = 0) text(C[, axespar], dimnames(C)[[1]]) } } cat("_______________________________________________________________\n" ) if(cor == F) { repeat { cat("graphique pour les axes (o/n) ?\n") plta <- scan(quiet=T,"", character(), 1) if((length(plta) == 0) | (plta == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltah <- scan(quiet=T,"", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltav <- scan(quiet=T,"", numeric(), 1) par(mfrow = c(1, 1), pty = "s") axespar <- c(pltah, pltav) plot(A[, axespar], xlab = paste("a", axespar[ 1], " ", round(valp[axespar[1]], digits = 4 ), "(", round(valp[axespar[1]]/inertiaACPVI * 100, digits = 2), "%)"), ylab = paste("a", axespar[2], " ", round(valp[axespar[2]], digits = 4), "(", round(valp[axespar[2]]/ inertiaACPVI * 100, digits = 2), "%)"), type = "n") abline(h = 0) abline(v = 0) text(A[, axespar], dimnames(A)[[1]]) } } } else { repeat { cat("cercle des correlations pour les variables de Y (o/n) ?\n" ) plta <- scan(quiet=T,"", character(), 1) if((length(plta) == 0) | (plta == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltah <- scan(quiet=T,"", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltav <- scan(quiet=T,"", numeric(), 1) par(mfrow = c(1, 1), pty = "s") axespar <- c(pltah, pltav) theta <- seq(0, 20, 0.05) x <- cos(theta) y <- sin(theta) plot(x, y, type = "l", xlab = paste("a", axespar[1], " ", round(valp[axespar[1]], digits = 4), "(", round(valp[axespar[1]]/ inertiaACPVI * 100, digits = 2), "%)"), ylab = paste("a", axespar[2], " ", round( valp[axespar[2]], digits = 4), "(", round( valp[axespar[2]]/inertiaACPVI * 100, digits = 2), "%)"), type = "n") abline(h = 0) abline(v = 0) text(A[, axespar], dimnames(A)[[1]]) } } repeat { cat("cercle des correlations pour les variables de X (o/n) ?\n" ) plta <- scan(quiet=T,"", character(), 1) if((length(plta) == 0) | (plta == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltah <- scan(quiet=T,"", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltav <- scan(quiet=T,"", numeric(), 1) par(mfrow = c(1, 1), pty = "s") axespar <- c(pltah, pltav) theta <- seq(0, 20, 0.05) x <- cos(theta) y <- sin(theta) plot(x, y, type = "l", xlab = paste("a", axespar[1], " ", round(valp[axespar[1]], digits = 4), "(", round(valp[axespar[1]]/ inertiaACPVI * 100, digits = 2), "%)"), ylab = paste("a", axespar[2], " ", round( valp[axespar[2]], digits = 4), "(", round( valp[axespar[2]]/inertiaACPVI * 100, digits = 2), "%)"), type = "n") abline(h = 0) abline(v = 0) text(t(ZA[axespar, ]/sqrt(valp[axespar])), dimnames(ZA)[[2]]) } } } } return(list(nomfichX=nomfichX,Xini=Xini,X=X,QX=QX,nomfichY=nomfichY,Yini=Yini,Y=Y,QY=QY,D=D,inertiaACPVI=inertiaACPVI,k=k, valp=valp,A=A,C=C,ZA=ZA)) } ###################################################### acpxqd<-function(X, Q = 1, D = 1, centrer = T, cor = T,k=3, impres = T, graph = T, Xl = 1, Xc = 1, aideus = 1, aideva = 1,colpar=3,bgpar="wheat",askpar=F,cexpar=0.7) { # # DVS d'un triplet (X,Q,D) # # Entrees # # X matrice des variables # Q metrique des u.s. : si Q=1 alors Q=Idn, si Q est un vecteur alors Q=diag(vecteur), sinon Q=matrice # D metrique des variables, si D=1 alors 1/nIdn, si D est un vecteur alors D=diag(vecteur), sinon D=matrice # centrer=F on ne centre pas X, sinon centrer=T # cor=F on ne reduit pas, centrage et reduction si cor=T # k nombre de composantes retenues, par defaut 3 # impres impression des resultats si T, si F rien # graph trace des graphiques si T, si F pas de trace # Xl matrice eventuelle des u.s. supplememtaires # Xc matrice eventuelle des variables supplememtaires # aideus=1 aides a l'interpretation pour les u.s., CTR et COS, pour les k premiers axes, sinon aideus=0 # aideva=1 aides a l'interpretation pour les variables, CTR et COS, pour les k premiers axes, sinon aideva=0 # bgpar chaine de caractères, couleur du background, par exple, "bisque1", "whitesmoke"...utiliser colors() pour choisir # Sorties # # nomfichX nom du fichier X # Xini tableau initial # X tableau apres eventuel centrage et/ou reduction # Q metrique utilisee pour les us # D metrique utilisee pour les variables # inertiaX inertie du triplet # k rang de la DSD retenue # valp vecteur de toutes les valeurs propres # Fa matrice des k premiers facteurs principaux de la DSD # A matrice des k premiers axes principaux de la DSD # C matrice des k premieres composantes principales de la DSD # CRTus et COSus contributions absolues et relatives des u.s. pour les k premieres composantes # CTRva et COSva contributions absolues et relatives des variables pour les k premiers axes # As matrice des k premiers axes pour les eventuelles variables supplementaires # Cs matrice des k premieres composantes pour les eventuelles u.s. supplementaires # # Fonctions exterieures appelees # # Dcentred # # R.Sabatier revu par JF Durand # MAJ 03/04/97 20/01/2002 # lecture des donnees nomfichX <- deparse(substitute(X)) X <- as.matrix(X) # if(length(c(which.inf(X), which.na(X)))) # stop("valeurs manquantes ou infinies dans X") if(impres){ oldpar<-par(no.readonly = TRUE) par(ask=F,bg=bgpar,mfrow=c(1,1),bty="n",xaxt="n",yaxt="n") plot(0,0,type="n",bg=par("bg"),xlab="",ylab="",xlim=c(-4,4),ylim=c(-4,4)) text(0,2,"Anlyse en Composantes",cex=1.8,col="red") text(0,0,"Principales du triplet",cex=1.8,col="red") text(0,-2,"( X , Q , D )",cex=1.8,col="red") par(oldpar) par(ask=askpar,bg=bgpar,xaxt="s",yaxt="s") } n <- nrow(X) p <- ncol(X) if(is.null(dimnames(X))) dimnames(X) <- list(format(1:n), paste("v", 1:p, sep = "")) if(length(dimnames(X)[[1]]) == 0) dimnames(X)[[1]] <- format(1:n) if(length(dimnames(X)[[2]]) == 0) dimnames(X)[[2]] <- paste("v", 1:p, sep = "") Xini <- as.matrix(X) if(p < 2) return() if(impres == F) graph <- F if(impres == T) { if(!exists(".Device", frame = 0)) { cat("Initialisez le graphique !!!\n") return() } } if(centrer == F) cor <- F if(cor == T) centrer <- T p2 <- p * p n2 <- n * n As <- NULL Cs <- NULL CTRus <- NULL COSus <- NULL CTRva <- NULL COSva <- NULL # calcul de la metrique Q if(length(Q) == 1) Q <- diag(p) if(length(Q) == p) Q <- diag(Q, nrow = p) if(length(Q) == p2) Q <- as.matrix(Q) # calcul de la metrique D if(length(D) == 1) D <- diag(rep(1/n, n), nrow = n) if(length(D) == n) D <- diag(D, nrow = n) if(length(D) == n2) D <- as.matrix(D) # verifications u.s. et variables supplementaires if(length(Xl) != 1) { isup <- 1 Xl <- as.matrix(Xl) if(ncol(Xl) != ncol(X)) return() if(is.null(dimnames(Xl))) dimnames(Xl) <- list(paste("is", 1:n, sep = ""), paste( "v", 1:p, sep = "")) if(length(dimnames(Xl)[[1]]) == 0) dimnames(Xl)[[1]] <- paste("is", 1:n, sep = "") if(length(dimnames(Xl)[[2]]) == 0) dimnames(Xl)[[2]] <- paste("v", 1:p, sep = "") } else isup <- 0 if(length(Xc) != 1) { vsup <- 1 Xc <- as.matrix(Xc) if(nrow(Xc) != nrow(X)) return() if(is.null(dimnames(Xc))) dimnames(Xc) <- list(paste("i", 1:n, sep = ""), paste( "vs", 1:p, sep = "")) if(length(dimnames(Xc)[[1]]) == 0) dimnames(Xc)[[1]] <- paste("i", 1:n, sep = "") if(length(dimnames(Xc)[[2]]) == 0) dimnames(Xc)[[2]] <- paste("vs", 1:p, sep = "") } else vsup <- 0 if(centrer == T) { centrage <- Dcentred(X, D = D) X <- centrage$Xc meanX <- matrix(centrage$moy) varX <- matrix(centrage$var) dimnames(meanX) <- list(dimnames(X)[[2]], "moy") dimnames(varX) <- list(dimnames(X)[[2]], "var") if(isup == 1) Xl <- sweep(Xl, 2, meanX) if(vsup == 1) { centragevs <- Dcentred(Xc, D = D) Xc <- centragevs$Xc } } if(cor == T) { X <- centrage$Xcr if(isup == 1) Xl <- sweep(Xl, 2, sqrt(varX), FUN = "/") if(vsup == 1) Xc <- centragevs$Xcr } if(impres) { cat(" - D.V.S. d'un triplet (X,Q,D) -\n") cat(" -------------------------------\n") if(isup == 1) cat(" avec u.s. supplementaires\n") if(vsup == 1) cat(" avec variables supplementaires\n") cat("__________________________________________________________________\n" ) if(cor == T) { cat(" donnees centrees et reduites\n") cat(" moyenne des variables de X\n") print(t(meanX)) cat(" variance des variables de X\n") print(t(varX)) } if(centrer == T && cor == F) { cat(" donnees centrees\n") cat(" moyenne des variables de X\n") print(t(meanX)) cat(" variance des variables de X\n") print(t(varX)) } cat("__________________________________________________________________\n" ) } # initialisation valp <- NULL # vecteur des valeurs propres A <- NULL # matrice des facteurs principaux C <- NULL # matrice des composantes principales # decomposition de Choleski de Q Qhalf <- chol(round(Q, 6)) # Qhalf <- chol(round(Q, 6), pivot = T) # Qhalf <- Qhalf[, order(attr(Qhalf, "pivot"))] Y <- X %*% t(Qhalf) # diagonalisation diago <- t(Y) %*% D %*% Y eg <- eigen(diago, symmetric = T) valp <- eg$values if(all(eg$vectors[,1]<=0)) eg$vectors[,1]=-eg$vectors[,1] for(i in 1:ncol(diago)) if(valp[i] < 0) valp[i] <- 0 inertiaX <- sum(valp) valp2 <- sum(valp^2) # impression de l'histogramme des valeurs propres if(impres) { par(bg=bgpar) cat(paste("Inertie totale de (X,Q,D) =", format(inertiaX), "\n" )) # choix du nombre d'axes si non donne en argument cat("__________________________________________________________________\n" ) valtab <- matrix(0, p, 4) dimnames(valtab) <- list(format(1:p), c("val.pro.", "% inert.", "% cumul.", " RV")) valtab[, 1] <- round(valp, digits = 4) valtab[, 2] <- round(valp/inertiaX * 100, digits = 2) valtab[, 3] <- round(cumsum(valp/inertiaX * 100),2) for(i in 1:p) { #valtab[i, 3] <- sum(valtab[1:i, 2]) valtab[i, 4] <- sqrt(sum(valp[1:i]^2)/valp2) } print(valtab) cat("barplots des valeurs propres (o ou n) ?\n") plth <- scan(quiet=T,"", character(), 1) if(length(plth) == 0) plth <- "n" if(plth == "o" || plth == "O") { par(mfrow = c(1, 2),pty="m") barplot(valp, space = 2, names = format(1:p),bg=par("bg"),cex.names =cexpar,cex.axis=cexpar,col=palette()[1:length(valp)+1]) title("valeurs propres",cex.main=cexpar+0.4) barplot(valtab[,3], space = 2, names = format(1:p),bg=par("bg"),cex.names =cexpar,cex.axis=cexpar,col=palette()[1:length(valp)+1]) abline(h=100,lty=3,col="red") title("% inertie cumulée",cex.main=cexpar+0.4) } repeat { cat("Combien de composantes voulez-vous ? (<=", p, ")") k <- scan(quiet=T,"", n = 1) if(k != 0) break } } # calcul des composantes principales Cent <- Y %*% eg$vectors C <- Cent[, 1:k] if(k == 1) C <- as.matrix(C) dimnames(C) <- list(dimnames(X)[[1]], paste("c", format(1:k), sep = "") ) # calcul des axes principaux Aent <- matrix(NA, nrow = ncol(X), ncol = p) if(length(Q) != 1) fac <- solve(Qhalf) %*% eg$vectors else fac <- eg$vectors for(i in 1:p) Aent[, i] <- (fac[, i]) * sqrt(valp[i]) A <- Aent[, 1:k] dimnames(A) <- list(dimnames(X)[[2]], paste("a", format(1:k), sep = "") ) # calcul des facteurs principaux Fa <- matrix(NA, nrow = ncol(X), ncol = k) Fa <- Q %*% A dimnames(Fa) <- list(dimnames(X)[[2]], paste("f", format(1:k), sep = "" )) # calcul des coordonnees u.s. et variables supplementaires if(isup == 1) { Cs <- Xl %*% Q %*% A for(i in 1:k) Cs[, i] <- (Cs[, i])/sqrt(valp[i]) dimnames(Cs) <- list(dimnames(Xl)[[1]], paste("c", format(1:k), sep = "")) } if(vsup == 1) { As <- t(Xc) %*% D %*% C for(i in 1:k) As[, i] <- (As[, i])/sqrt(valp[i]) dimnames(As) <- list(dimnames(Xc)[[2]], paste("a", format(1:k), sep = "")) } # calcul et impression des aides a l'interpretation if(aideus == 1) { CTRus <- matrix(NA, nrow = n, ncol = k) COSus <- CTRus for(i in 1:n) for(j in 1:k) CTRus[i, j] <- round((D[i, i] * C[i, j]^2)/valp[ j] * 10000, digits = 0) dimnames(CTRus) <- list(dimnames(X)[[1]], paste("CTA", format(1: k), sep = "")) for(i in 1:n) { ss <- 0 for(j in 1:p) ss <- ss + Cent[i, j]^2 for(j in 1:k) COSus[i, j] <- round(Cent[i, j]^2/ss * 10000, digits = 0) } dimnames(COSus) <- list(dimnames(X)[[1]], paste("COS", format(1: k), sep = "")) } if(aideva == 1) { CTRva <- matrix(NA, nrow = p, ncol = k) COSva <- CTRva for(i in 1:p) for(j in 1:k) CTRva[i, j] <- round((Q[i, i] * A[i, j]^2)/valp[ j] * 10000, digits = 0) dimnames(CTRva) <- list(dimnames(X)[[2]], paste("CTA", format(1: k), sep = "")) for(i in 1:p) { ss <- 0 for(j in 1:p) ss <- ss + Aent[i, j]^2 for(j in 1:k) COSva[i, j] <- round(Aent[i, j]^2/ss * 10000, digits = 0) } dimnames(COSva) <- list(dimnames(X)[[2]], paste("COS", format(1: k), sep = "")) } if(impres) { if(aideus == 1) { cat("_______________________________________________________________\n" ) cat("aides a l'interpretation pour les u.s. (o/n) ?\n") plta <- scan(quiet=T,"", character(), 1) if(length(plta) == 0) plta <- "n" if(plta == "o" || plta == "O") { cat("Contributions absolues des", n, "u.s. pour les", k, "premieres composantes\n" ) print(CTRus) cat("Contributions relative des", n, "u.s. pour les", k, "premieres composantes\n" ) print(COSus) } } if(aideva == 1) { cat("_______________________________________________________________\n" ) cat("aides a l'interpretation pour les variables (o/n) ?\n" ) plta <- scan(quiet=T,"", character(), 1) if(length(plta) == 0) plta <- "n" if(plta == "o" || plta == "O") { cat("Contributions absolues des", p, "variables pour les", k, "premiers axes\n") print(CTRva) cat("Contributions relative des", p, "variables pour les", k, "premiers axes\n") print(COSva) } } } # trace des eventuels graphiques if(graph == T) { cat("_______________________________________________________________\n" ) repeat { cat("graphique pour les u.s. (o/n) ?\n") pltc <- scan(quiet=T,"", character(), 1) if(length(pltc) == 0) pltc<-"n" if((pltc == "N") | (pltc == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltch <- scan(quiet=T,"", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltcv <- scan(quiet=T,"", numeric(), 1) par(mfrow = c(1, 1), pty = "s") axespar <- c(pltch, pltcv) if(isup == 1) Ctot <- rbind(C, Cs) else Ctot <- C plot(Ctot[, axespar], xlab = paste("Dim", axespar[ 1], "(", round(valp[axespar[1]]/inertiaX * 100, digits = 2), "%)"), ylab = paste("Dim", axespar[ 2],"(", round(valp[axespar[2]]/inertiaX * 100, digits = 2), "%)"), type = "n",bg=par("bg"),cex=cexpar) abline(h = 0) abline(v = 0) text(Ctot[, axespar], dimnames(Ctot)[[1]],col=colpar+1,cex=cexpar) cat("Une carte des individus actifs coloriés selon les modalités d'une variable qualitative? (o/n)\n") plto <- scan(quiet=T,"", character(), 1) if( (length(plto)!=0)&&((plto=="O")|(plto=="o"))) { repeat{ cat("entrer le nom du vecteur-variable ou celui de la matrice contenant la colonne-variable\n") nom <- scan(quiet=T,"", character(), 1) matricenom<-as.matrix(get(nom,pos=1)) if(nrow(matricenom)==nrow(Xini))break else cat("misfit du nombre de lignes\n") } if(is.vector(get(nom,pos=1))) {dimnames(matricenom)<-list(dimnames(Xini)[[1]],c(nom)) numer<-1 variab<-get(nom,pos=1) } if((is.matrix(get(nom,pos=1)))|(is.data.frame(get(nom,pos=1)))) { cat(paste(dimnames(matricenom)[[2]],"(",1:ncol(matricenom),"),",sep=""),"\n") cat("entrer le numero de la colonne de valeurs entieres\n") numer <- scan(quiet=T,"",numeric(), 1) variab<-matricenom[,numer] cat("nom de la variable : ",dimnames(matricenom)[[2]][numer],"\n") } cat("click to locate the top left corner of the legend\n") for(i in min(variab):max(variab)) text(C[variab==i,axespar],dimnames(C)[[1]][variab==i],col=i+1,cex=cexpar) legend(locator(1),paste(dimnames(matricenom)[[2]][numer],min(variab):max(variab),sep=""), fill=(min(variab):max(variab))+1,cex=cexpar) } } } cat("_______________________________________________________________\n" ) repeat { if(cor == F) { cat("graphique pour les variables (o/n) ?\n") plta <- scan(quiet=T,"", character(), 1) if(length(plta) == 0) plta<-"n" if((plta == "N") | (plta == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltah <- scan(quiet=T,"", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltav <- scan(quiet=T,"", numeric(), 1) par(mfrow = c(1, 1), pty = "s") axespar <- c(pltah, pltav) if(vsup == 1) Atot <- rbind(A, As) else Atot <- A plot(rbind(Atot[, axespar],c(0,0)), xlab = paste("Dim", axespar[1]," (", round(valp[axespar[1]]/ inertiaX * 100, digits = 2), "%)"), ylab = paste("Dim", axespar[2]," (", round(valp[ axespar[2]]/inertiaX * 100, digits = 2), "%)"), type = "n",bg=par("bg"),cex=cexpar) arrows(rep(0,p),rep(0,p),0.92*Atot[,axespar[1]],0.92*Atot[,axespar[2]],col=colpar+3,angle=15,length=0.15) abline(v=0) abline(h=0) text(Atot[, axespar], dimnames(Atot)[[1]],col=colpar+3,cex=cexpar) if(vsup==1){ arrows(rep(0,nrow(As)),rep(0,nrow(As)),0.92*As[,axespar[1]],0.92*As[,axespar[2]],col=colpar+2,angle=15,length=0.15) text(As[, axespar,drop=F], dimnames(Xc)[[2]],col=colpar+2) } } } else { cat("cercle des correlations (o/n) ?\n") plta <- scan(quiet=T,"", character(), 1) if(length(plta) == 0) plta<-"n" if((plta == "N") | (plta == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltah <- scan(quiet=T,"", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltav <- scan(quiet=T,"", numeric(), 1) par(mfrow = c(1, 1), pty = "s") axespar <- c(pltah, pltav) theta <- seq(0, 20, 0.05) x <- cos(theta) y <- sin(theta) if(vsup == 1) Atot <- rbind(A, As) else Atot <- A plot(x, y, type = "l", xlab = paste("Dim", axespar[1], " (", round(valp[axespar[1]]/ inertiaX * 100, digits = 2), "%)"), ylab = paste("Dim", axespar[2]," (", round(valp[ axespar[2]]/inertiaX * 100, digits = 2), "%)"),cex=cexpar) abline(h = 0) abline(v = 0) arrows(rep(0,p),rep(0,p),0.92*Atot[, axespar[1]],0.92*Atot[, axespar[2]],col=colpar+3,angle=15,length=0.15) text(Atot[, axespar], dimnames(Atot)[[1]],col=colpar+3,cex=cexpar) if(vsup==1){ arrows(rep(0,nrow(As)),rep(0,nrow(As)),0.92*As[,axespar[1]],0.92*As[,axespar[2]],col=colpar+2,angle=15,length=0.15) text(As[, axespar,drop=F], dimnames(Xc)[[2]],col=colpar+2,cex=cexpar) } } } } # fin des graphiques } invisible(return(list(nomfichX=nomfichX,Xini=Xini,X=X,Q=Q,D=D,inertiaX=inertiaX,k=k,valp=valp,Fa=Fa,A=A,C=C,CTRus=CTRus, COSus=COSus,CTRva=CTRva,COSva=COSva,As=As,Cs=Cs))) } ############################################## afc<-function(N, AFCM=F,lignesup,colsup,DI = 1, DJ = 1, M = 1, k = 3, ns = F, impres = T, graph = T,colpar=1, bgpar="lightblue",cexpar=0.7,askpar=F) { # # afc et afc generalisee par rapport a un modele avec metriques lignes et colonnes quelconques # # afcm si AFCM=T, alors N est le resultat de MVcut(X) ou de Codisjc(X) donnant le codage disjonctif # complet du tabeau X. Exemple: afc(MVcut(X),AFCM=T) # # # Entrees # # N tableau de contingence ou resultat de MVcut ou de Codisjc si AFCM=T, # lignesup, matrice des lignes supplementaires # colsup, matrice des colonnes supplementaires # DI metrique des lignes, si absent c'est la marge # DJ metrique des colonnes, si absent c'est la marge # M modele (matrice meme dim que N) si c'est une afc par rapport a un modele, sinon M=1 # k nombre de composantes retenues, par defaut 3 # ns=T si afc non symetrique, si ns=F afc normale # impres impression des resultats si T, si F rien # graph trace des graphiques si T, si F pas de trace # bgpar couleur du background graphique, utiliser aussi "wheat" ou "lightgray" utiliser colors() # pour choisir # Sorties # # nomfichX nom du fichier X # Nini tableau initial # X tableau apres eventuelle transformation # DI metrique ligne # DJ metrique colonne # inertia inertie du triplet # k rang de la DSD retenue # valp vecteur de toutes les valeurs propres # A matrice des k premiers axes de la DSD # C matrice des k premieres composantes principales de la DSD # CTRus et COSus contributions absolues et relatives des u.s. pour les k premieres composantes # CTRva et COSva contributions absolues et relatives des variables pour les k premiers axes # # Fonctions exterieures appelees # # acpxqd, Dcentred # # R.Sabatier revu par JF Durand # MAJ 03/04/97 25/11/2003 # # lecture des donnees nomfichN <- deparse(substitute(N)) oldpar<-par(no.readonly = TRUE) par(ask=F,bg=bgpar,mfrow=c(1,1),bty="n",xaxt="n",yaxt="n") plot(0,0,type="n",bg=par("bg"),xlab="",ylab="",xlim=c(-4,4),ylim=c(-4,4)) text(0,2,"Analyse des Correspondances",cex=1.7,col="red") if(AFCM)text(0,0,"Multiples",cex=1.7,col="red") else text(0,0,"Simples",cex=1.7,col="red") par(oldpar) par(ask=askpar,bg=bgpar,xaxt="s",yaxt="s") if(AFCM){ nbmod<-N$nbmod labels<-N$labels N<-N$U } N <- as.matrix(N) M <- as.matrix(M) # if(length(c(which.inf(N), which.na(N)))) # stop("valeurs manquantes ou infinies dans N") # if(length(c(which.inf(M), which.na(M)))) # stop("valeurs manquantes ou infinies dans M") Nini <- N I <- nrow(N) J <- ncol(N) IJ <- I * J if(I <= 2) return() if(J <= 2) return() if(length(M) != 1 && ns == T) stop("conflits de parametres entre M et ns") if(impres == F) graph <- F n <- sum(N) P <- N/n I1 <- matrix(1, nrow = I, ncol = 1) J1 <- matrix(1, nrow = J, ncol = 1) # metriques si afc non symetrique if(ns == T) { DJ <- diag(J) DI <- diag(P %*% J1, nrow = I) } else { # calcul de la metrique DI if(length(DI) == 1) DI <- diag(as.vector(P %*% J1)) if(length(DI) == I) DI <- diag(as.vector(DI)) # calcul de la metrique DJ if(length(DJ) == 1) DJ <- diag(as.vector(t(I1) %*% P)) if(length(DJ) == J) DJ <- diag(as.vector(DJ)) } # les elements supplementaires doivent etre des frequences if(!missing(lignesup)) {lignesup<-as.matrix(lignesup) lignesup<-sweep(lignesup,1,apply(lignesup,1,sum),FUN="/")} if(!missing(colsup)) {colsup<-as.matrix(colsup) colsup<-sweep(colsup,2,apply(colsup,2,sum),FUN="/")} if(impres) { cat("_______________________________________________________________\n" ) if(ns == F) cat(" - afc de ", nomfichN, " -\n") if(ns == T) cat(" - afc non symetrique de ", nomfichN, " -\n") if(length(M) == IJ) cat(" par rapport a un modele\n") } # calcul de X, mis dans N if(length(M) == 1) { if(ns == F) for(i in 1:I) for(j in 1:J) N[i, j] <- P[i, j]/(DJ[j, j] * DI[i, i]) - 1 if(ns == T) for(i in 1:I) for(j in 1:J) N[i, j] <- (P[i, j]/DI[i, i]) - diag(t(I1) %*% P, nrow = J)[j, j] } if(length(M) == IJ) for(i in 1:I) for(j in 1:J) N[i, j] <- (P[i, j] - M[i, j]/n)/(DJ[j, j] * DI[ i, i]) res <- acpxqd(N, Q = DJ, D = DI, centrer = F, k = k, impres = F, graph = F) inertia <- res$inertiaX chi2 <- n * inertia ddl <- (I - 1) * (J - 1) valp <- res$valp C <- res$C k <- res$k A <- res$A if(!missing(colsup))coordcolsup<-t(colsup)%*%C[,1:k]%*%solve(sqrt(diag(valp[1:k]))) if(!missing(lignesup)) coordlignesup<-lignesup%*%A[,1:k]%*%solve(sqrt(diag(valp[1:k]))) CTRus <- res$CTRus COSus <- res$COSus CTRva <- res$CTRva COSva <- res$COSva if(impres == T) { cat("_______________________________________________________________\n" ) cat(paste("Inertie totale =", format(inertia), "\n")) cat(paste("nb de lignes =", I, "\n")) cat(paste("nb de colonnes =", J, "\n")) cat(paste("effectif total =", n, "\n")) if((ns == F)&(AFCM==F)) cat(paste("d2 d'indépendance=", format(chi2 ),",d.d.l.=",ddl, ",Chi2 crit.=", format( qchisq(0.95, ddl)),"(0.05)",sep=""), "\n") cat("_______________________________________________________________\n" ) pp <- length(valp) valtab <- matrix(0, pp, 3) dimnames(valtab) <- list(format(1:pp), c("val.pro.", "% inert.", "% cumul.")) valtab[, 2] <- valp/inertia * 100 for(i in 1:pp) valtab[i, 3] <- sum(valtab[1:i, 2]) valtab[, 1] <-round(res$valp,5) valtab[, 2:3] <-round(valtab[,2:3],2) print(valtab) if(graph){ cat("diagrammes des valeurs propres (o ou n) ?") plth <- scan(quiet=T,"", character(), 1) if(length(plth) == 0) plth <- "n" if(plth == "o" || plth == "O") { par(bg=bgpar,mfrow=c(1,2),pty="m") palette(rainbow(12)) barplot(valp, space = 2, names = format(1:J), col = 1:max(I,J), cex.names =cexpar,cex.axis=cexpar) title("valeurs propres",cex.main=cexpar+0.4) barplot(valtab[,3], space = 2, names = format(1:J), col = 1:max(I,J), cex.names =cexpar,cex.axis=cexpar) title("% inertie cumulée",cex.main=cexpar+0.4) abline(h=100,lty=3,col="red") # palette(rainbow(6)) } repeat { cat("Combien de composantes voulez-vous ? (<=",min(I,J), ")") k <- scan(quiet=T,"", n = 1) if(k != 0) break } } } # calcul et impression des aides a l'interpretation if((impres == T) | length(D) == I) { cat("_______________________________________________________________\n" ) cat("aides a l'interpretation pour les modalites lignes (o/n) ?\n" ) plta <- scan(quiet=T,"", character(), 1) if(length(plta) == 0) plta <- "n" if(plta == "o" || plth == "O") { cat("Contributions absolues des", I, "modalites pour les", k, "premieres composantes\n") print(CTRus) cat("\nContributions relative des", I, "modalites pour les", k, "premieres composantes\n") print(COSus) } } if((impres == T) | length(D) == J) { cat("_______________________________________________________________\n" ) cat("aides a l'interpretation pour les modalites colonnes (o/n) ?\n" ) plta <- scan(quiet=T,"", character(), 1) if(length(plta) == 0) plta <- "n" if(plta == "o" || plth == "O") { cat("Contributions absolues des", J, "modalites pour les", k, "premiers axes\n") print(CTRva) cat("\nContributions relative des", J, "modalites pour les", k, "premiers axes\n") print(COSva) } } # trace des eventuels graphiques if(graph == T) { cat("_______________________________________________________________\n" ) repeat { cat("graphique pour les composantes (lignes) (o/n) ?\n" ) pltc <- scan(quiet=T,"", character(), 1) if(length(pltc) == 0)pltc<-"n" if((pltc == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltch <- scan(quiet=T,"", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltcv <- scan(quiet=T,"", numeric(), 1) par(mfrow = c(1, 1), pty = "s") axespar <- c(pltch, pltcv) if(missing(lignesup)) plot(C[, axespar], xlab = paste("c", axespar[1], " ", round(valp[axespar[1]], digits = 4), "(", round(valp[axespar[1]]/inertia * 100, digits = 2), "%)"), ylab = paste("c", axespar[2], " ", round(valp[axespar[2]], digits = 4), "(", round(valp[axespar[2]]/inertia * 100, digits = 2), "%)"), type = "n",cex=cexpar) else plot(rbind(C[, axespar],coordlignesup[,axespar]), xlab = paste("c", axespar[1], " ", round(valp[axespar[1]], digits = 4), "(", round(valp[axespar[1]]/inertia * 100, digits = 2), "%)"), ylab = paste("c", axespar[2], " ", round(valp[axespar[2]], digits = 4), "(", round(valp[axespar[2]]/inertia * 100, digits = 2), "%)"), type = "n",cex=cexpar) abline(h = 0) abline(v = 0) text(C[, axespar], dimnames(C)[[1]], col = colpar+4,cex=cexpar) if(!missing(lignesup)) text(coordlignesup[,axespar,drop=F],dimnames(lignesup)[[1]],col=colpar+3,cex=cexpar) if(AFCM) { cat("Une carte des individus colories selon les modalites d'une variable qualitative? (o/n)\n") plto <- scan(quiet=T,"", character(), 1) if( (length(plto)!=0)&&((plto=="O")|(plto=="o"))) { repeat{ cat("entrer le nom du vecteur-variable ou celui de la matrice contenant la colonne-variable\n") nom <- scan(quiet=T,"", character(), 1) matricenom<-as.matrix(get(nom,pos=1)) if(nrow(matricenom)==nrow(C))break else cat("misfit du nombre de lignes\n") } if(is.vector(get(nom,pos=1))) {dimnames(matricenom)<-list(dimnames(C)[[1]],c(nom)) numer<-1 variab<-get(nom,pos=1) } if((is.matrix(get(nom,pos=1)))|(is.data.frame(get(nom,pos=1)))) { cat(paste(dimnames(matricenom)[[2]],"(",1:ncol(matricenom),"),",sep=""),"\n") cat("entrer le numero de la colonne de valeurs entieres\n") numer <- scan(quiet=T,"",numeric(), 1) variab<-matricenom[,numer] cat("nom de la variable : ",dimnames(matricenom)[[2]][numer],"\n") } cat("click to locate the top left corner of the legend\n") for(i in min(variab):max(variab)) text(C[variab==i,axespar],dimnames(matricenom)[[1]][variab==i],col=i+1,cex=cexpar) legend(locator(1),paste(dimnames(matricenom)[[2]][numer],min(variab):max(variab),sep=""), fill=(min(variab):max(variab))+1,cex=cexpar) } }#fin ifAFCM } } cat("_______________________________________________________________\n" ) repeat { cat("graphique pour les axes (colonnes) (o/n) ?\n") plta <- scan(quiet=T,"", character(), 1) if(length(plta) == 0)plta<-"n" if((plta == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltah <- scan(quiet=T,"", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltav <- scan(quiet=T,"", numeric(), 1) par(mfrow = c(1, 1), pty = "s") axespar <- c(pltah, pltav) if(missing(colsup)) plot(A[, axespar], xlab = paste("a", axespar[1], " ", round(valp[axespar[1]], digits = 4), "(", round(valp[axespar[1]]/inertia * 100, digits = 2), "%)"), ylab = paste("a", axespar[2], " ", round(valp[axespar[2]], digits = 4), "(", round(valp[axespar[2]]/inertia * 100, digits = 2), "%)"), type = "n",cex=cexpar) else plot(rbind(A[, axespar],coordcolsup[,axespar]), xlab = paste("a", axespar[1], " ", round(valp[axespar[1]], digits = 4), "(", round(valp[axespar[1]]/inertia * 100, digits = 2), "%)"), ylab = paste("a", axespar[2], " ", round(valp[axespar[2]], digits = 4), "(", round(valp[axespar[2]]/inertia * 100, digits = 2), "%)"), type = "n",cex=cexpar) abline(h = 0) abline(v = 0) text(A[, axespar], dimnames(A)[[1]], col = colpar,cex=cexpar) if(!missing(colsup)) text(coordcolsup[,axespar,drop=F],dimnames(colsup)[[2]],col=colpar+5,cex=cexpar) if(sum((diag(DI)-(1/I))^2)<1e-16){# cas AFCM repeat {#repeat cat("relier des modalités d'une variable (o/n) ?\n") plta <- scan(quiet=T,"", character(), 1) if(length(plta) == 0)plta<-"n" if((plta == "n")) break else { cat("numero de la variable (<=",length(nbmod),"), et de la couleur\n") num<-scan(quiet=T,"",numeric(),2) if(num[1]==1)cumul<-0 else cumul<-sum(nbmod[1:(num[1]-1)]) lines(A[(cumul+1):(cumul+nbmod[num[1]]),axespar],col=num[2]) text(A[(cumul+1):(cumul+nbmod[num[1]]),axespar],dimnames(A)[[1]][(cumul+1):(cumul+nbmod[num[1]])],col=num[2],cex=cexpar) } }#repeat }# } } cat("_______________________________________________________________\n" ) repeat { cat("representations simultanees (o/n) ?\n") plta <- scan(quiet=T,"", character(), 1) if(length(plta) == 0)plta<-"n" if((plta == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltah <- scan(quiet=T,"", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltav <- scan(quiet=T,"", numeric(), 1) par(mfrow = c(1, 1), pty = "s") axespar <- c(pltah, pltav) Z <- matrix(0, nrow = I + J, ncol = 2) Z[, 1] <- c(C[, axespar[1]], A[, axespar[1]]) Z[, 2] <- c(C[, axespar[2]], A[, axespar[2]]) dimnames(Z) <- list(c(dimnames(C)[[1]], dimnames(A)[[1]]), c("1", "2")) if(missing(colsup)& missing(lignesup)) Ztot<-Z else { if((!missing(colsup))&(!missing(lignesup))) Ztot<-rbind(Z,coordcolsup[,c(axespar[1],axespar[2])],coordlignesup[,c(axespar[1],axespar[2])]) if((!missing(colsup))&(missing(lignesup))) Ztot<-rbind(Z,coordcolsup[,c(axespar[1],axespar[2])]) if((missing(colsup))&(!missing(lignesup))) Ztot<-rbind(Z,coordlignesup[,c(axespar[1],axespar[2])]) } plot(Ztot, xlab = paste("axe",axespar[1], " ", round( valp[axespar[1]], digits = 4), "(", round( valp[axespar[1]]/inertia * 100, digits = 2), "%)"), ylab = paste("axe",axespar[2], " ", round( valp[axespar[2]], digits = 4), "(", round( valp[axespar[2]]/inertia * 100, digits = 2), "%)"), type = "n",cex=cexpar) abline(h = 0) abline(v = 0) text(Z[1:I, ], dimnames(Z)[[1]][1:I], col =colpar+ 4,cex=cexpar) text(Z[(I + 1):(I + J), ], dimnames(Z)[[1]][(I + 1):(I + J)], col =colpar,cex=cexpar) if(!missing(colsup)) text(coordcolsup[,c(axespar[1],axespar[2]),drop=F],dimnames(colsup)[[2]],col=colpar+5,cex=cexpar) if(!missing(lignesup)) text(coordlignesup[,c(axespar[1],axespar[2]),drop=F],dimnames(lignesup)[[1]],col=colpar+10,cex=cexpar) if(AFCM) { cat("Une carte des individus actifs colories selon les modalites d'une variable qualitative? (o/n)\n") plto <- scan(quiet=T,"", character(), 1) if( (length(plto)!=0)&&((plto=="O")|(plto=="o"))) { repeat{ cat("entrer le nom du vecteur-variable ou celui de la matrice contenant la colonne-variable\n") nom <- scan(quiet=T,"", character(), 1) matricenom<-as.matrix(get(nom,pos=1)) if(nrow(matricenom)==nrow(C))break else cat("misfit du nombre de lignes\n") } if(is.vector(get(nom,pos=1))) {dimnames(matricenom)<-list(dimnames(C)[[1]],c(nom)) numer<-1 variab<-get(nom,pos=1) } if((is.matrix(get(nom,pos=1)))|(is.data.frame(get(nom,pos=1)))) { cat(paste(dimnames(matricenom)[[2]],"(",1:ncol(matricenom),"),",sep=""),"\n") cat("entrer le numero de la colonne de valeurs entieres\n") numer <- scan(quiet=T,"",numeric(), 1) variab<-matricenom[,numer] cat("nom de la variable : ",dimnames(matricenom)[[2]][numer],"\n") } cat("click to locate the top left corner of the legend\n") for(i in min(variab):max(variab)) text(C[variab==i,axespar],dimnames(C)[[1]][variab==i],col=i+1,cex=cexpar) legend(locator(1),paste(dimnames(matricenom)[[2]][numer],min(variab):max(variab),sep=""), fill=(min(variab):max(variab))+1,cex=cexpar) } }#fin ifAFCM if(sum((diag(DI)-(1/I))^2)<1e-16){# repeat {#repeat cat("relier des modalités d'une variable (o/n) ?\n") plta <- scan(quiet=T,"", character(), 1) if(length(plta) == 0)plta<-"n" if((plta == "n")) break else {palette(rainbow(12)) cat("numero de la variable (<=",length(nbmod),"), et de la couleur\n") num<-scan(quiet=T,"",numeric(),2) if(num[1]==1)cumul<-0 else cumul<-sum(nbmod[1:(num[1]-1)]) lines(Z[(I+cumul+1):(I+cumul+nbmod[num[1]]),],col=num[2]) text(Z[(I+cumul+1):(I+cumul+nbmod[num[1]]),],dimnames(A)[[1]][(cumul+1):(cumul+nbmod[num[1]])],col=num[2],cex=cexpar) } palette(rainbow(6)) }#repeat }# else {##cas afc simple repeat {#repeat cat("relier des modalités d'une variable (o/n) ?\n") plta <- scan(quiet=T,"", character(), 1) if(length(plta) == 0)plta<-"n" if((plta == "n")) break else { cat("Rentrez 2 chiffres\n") cat("ligne (1) ou colonne (2) ou ligne sup (3) ou col sup(4) et le numero de la couleur\n") num<-scan(quiet=T,"",numeric(),2) if(num[1]==1){ lines(C[,axespar],col=num[2]) text(C[,axespar],dimnames(C)[[1]],col=num[2],cex=cexpar) } if(num[1]==2) { lines(A[,axespar],col=num[2]) text(A[,axespar],dimnames(A)[[1]],col=num[2],cex=cexpar)} if(num[1]==3) { lines(coordlignesup[,axespar],col=num[2]) text(coordlignesup[,axespar],dimnames(lignesup)[[1]],col=num[2],cex=cexpar)} if(num[1]==4) { lines(coordcolsup[,axespar],col=num[2]) text(coordcolsup[,axespar],dimnames(colsup)[[2]],col=num[2],cex=cexpar)} } }#repeat }## } } if(ns == F) { cat("_______________________________________________________________\n" ) repeat { cat("representations barycentriques (o/n) ?\n") plta <- scan(quiet=T,"", character(), 1) if(length(plta) == 0)plta<-"n" if((plta == "n")) break else { cat( "lignes au barycentre des colonnes (1) ou l'inverse (2) ?\n" ) bar <- scan(quiet=T,"", numeric(), 1) cat("axe horizontal (<=", k, ") ?\n") pltah <- scan(quiet=T,"", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltav <- scan(quiet=T,"", numeric(), 1) par(mfrow = c(1, 1), pty = "s") axespar <- c(pltah, pltav) if(bar == 1) { Z <- matrix(0, nrow = I + J, ncol = 2) Z[, 1] <- c(C[, axespar[1]], A[, axespar[1] ]/sqrt(valp[axespar[1]])) Z[, 2] <- c(C[, axespar[2]], A[, axespar[2] ]/sqrt(valp[axespar[2]])) dimnames(Z) <- list(c(dimnames(C)[[1]], dimnames(A)[[1]]), c("1", "2")) } if(bar == 2) { Z <- matrix(0, nrow = I + J, ncol = 2) Z[, 1] <- c(C[, axespar[1]]/sqrt(valp[ axespar[1]]), A[, axespar[1]]) Z[, 2] <- c(C[, axespar[2]]/sqrt(valp[ axespar[2]]), A[, axespar[2]]) dimnames(Z) <- list(c(dimnames(C)[[1]], dimnames(A)[[1]]), c("1", "2")) } plot(Z, xlab = paste("axe",axespar[1], " ", round( valp[axespar[1]], digits = 4), "(", round( valp[axespar[1]]/inertia * 100, digits = 2), "%)"), ylab = paste("axe",axespar[2], " ", round( valp[axespar[2]], digits = 4), "(", round( valp[axespar[2]]/inertia * 100, digits = 2), "%)"), type = "n",cex=cexpar) abline(h = 0) abline(v = 0) text(Z[1:I, ], dimnames(Z)[[1]][1:I], col = colpar+4,cex=cexpar) text(Z[(I + 1):(I + J), ], dimnames(Z)[[1]][( I + 1):(I + J)], col = colpar,cex=cexpar) if(AFCM) { cat("Une carte des individus actifs coloriés selon les modalités d'une variable qualitative? (o/n)\n") plto <- scan(quiet=T,"", character(), 1) if( (length(plto)!=0)&&((plto=="O")|(plto=="o"))) { repeat{ cat("entrer le nom du vecteur-variable ou celui de la matrice contenant la colonne-variable\n") nom <- scan(quiet=T,"", character(), 1) matricenom<-as.matrix(get(nom,pos=1)) if(nrow(matricenom)==nrow(C))break else cat("misfit du nombre de lignes\n") } if(is.vector(get(nom,pos=1))) {dimnames(matricenom)<-list(dimnames(C)[[1]],c(nom)) numer<-1 variab<-get(nom,pos=1) } if((is.matrix(get(nom,pos=1)))|(is.data.frame(get(nom,pos=1)))) { cat(paste(dimnames(matricenom)[[2]],"(",1:ncol(matricenom),"),",sep=""),"\n") cat("entrer le numero de la colonne de valeurs entieres\n") numer <- scan(quiet=T,"",numeric(), 1) variab<-matricenom[,numer] cat("nom de la variable : ",dimnames(matricenom)[[2]][numer],"\n") } cat("click to locate the top left corner of the legend\n") for(i in min(variab):max(variab)) text(Z[(1:I)[variab==i], ],dimnames(C)[[1]][variab==i],col=i+1,cex=cexpar) legend(locator(1),paste(dimnames(matricenom)[[2]][numer],min(variab):max(variab),sep=""), fill=(min(variab):max(variab))+1,cex=cexpar) } }#fin ifAFCM if(sum((diag(DI)-(1/I))^2)<1e-16){# repeat {#repeat cat("relier des modalités d'une variable (o/n) ?\n") plta <- scan(quiet=T,"", character(), 1) if(length(plta) == 0)plta<-"n" if((plta == "n")) break else {palette(rainbow(12)) cat("numero de la variable (<=",length(nbmod),"), et de la couleur\n") num<-scan(quiet=T,"",numeric(),2) if(num[1]==1)cumul<-0 else cumul<-sum(nbmod[1:(num[1]-1)]) lines(Z[(I+cumul+1):(I+cumul+nbmod[num[1]]),],col=num[2]) text(Z[(I+cumul+1):(I+cumul+nbmod[num[1]]),],dimnames(A)[[1]][(cumul+1):(cumul+nbmod[num[1]])],col=num[2],cex=cexpar) } palette(rainbow(6)) }#repeat }# } } } # fin des graphiques } palette("default") if((!missing(lignesup))& (!missing(colsup))) return(list(nomfichN=nomfichN,Nini=Nini,N=N,DI=DI,DJ=DJ,inertia=inertia,k=k,valp=valp,A=A,C=C,CTRus=CTRus,COSus=COSus, CTRva=CTRva,COSva=COSva,coordlignesup=coordlignesup,coordcolsup=coordcolsup)) if((!missing(colsup))& missing(lignesup))return(list(nomfichN=nomfichN,Nini=Nini,N=N,DI=DI,DJ=DJ,inertia=inertia,k=k,valp=valp,A=A,C=C,CTRus=CTRus,COSus=COSus, CTRva=CTRva,COSva=COSva,coordcolsup=coordcolsup)) if((!missing(lignesup))& missing(colsup))return(list(nomfichN=nomfichN,Nini=Nini,N=N,DI=DI,DJ=DJ,inertia=inertia,k=k,valp=valp,A=A,C=C,CTRus=CTRus,COSus=COSus, CTRva=CTRva,COSva=COSva,coordlignesup=coordlignesup)) if(missing(colsup)& missing(lignesup))return(list(nomfichN=nomfichN,Nini=Nini,N=N,DI=DI,DJ=DJ,inertia=inertia,k=k,valp=valp,A=A,C=C,CTRus=CTRus,COSus=COSus, CTRva=CTRva,COSva=COSva)) } ##################################################### posmultclas<-function(X, D = 1, dissim = T, k = 3, impres = T, graph = T, aideus = 1) { # # Positionement Multidimensionnel Classique # # Entrees # # X matrice des similarites ou dissimilarites # D metrique des variables, si D=1 alors 1/nIdn, si D est un vecteur alors D=diag(vecteur), sinon D=matrice # dissim=T la matrice est de dissimilarite, sinon dissim=F c'est une dissimilarite # k nombre de composantes retenues, par defaut 3 # impres impression des resultats si T, si F rien # graph trace des graphiques si T, si F pas de trace # aideus=1 aides a l'interpretation pour les u.s., CTR et COS, pour les k premiers axes, sinon aideus=0 # # Sorties # # nomfichX nom du fichier X # tableau initial # X tableau apres eventuel centrage et/ou reduction # D metrique des u.s. # W matrice de Torgerson # inertiaX inertie du triplet # k rang de la DSD retenue # valp vecteur de toute les valeurs propres # C matrice des k premieres composantes principales de la DSD # CRTus et COSus contributions absolues et relatives des u.s. pour les k premieres composantes # # lecture des donnees nomfichX <- deparse(substitute(X)) if(!is.matrix(X)) stop("X n'est pas une matrice") X <- as.matrix(X) if(length(c(which.inf(X), which.na(X)))) stop("valeurs manquantes ou infinies dans X") n <- nrow(X) if(is.null(dimnames(X))) dimnames(X) <- list(paste("i", 1:n, sep = ""), paste("v", 1:n, sep = "")) if(length(dimnames(X)[[1]]) == 0) dimnames(X)[[1]] <- paste("i", 1:n, sep = "") if(length(dimnames(X)[[2]]) == 0) dimnames(X)[[2]] <- paste("i", 1:n, sep = "") Xini <- as.matrix(X) if(n < 2) return() if(impres == F) graph <- F if(impres == T) { if(!exists(".Device", frame = 0)) { cat("Initialisez le graphique !!!\n") return() } } n2 <- n * n Cs <- NULL CTRus <- NULL COSus <- NULL if(length(D) > n) aideus <- 0 # calcul de la metrique D if(length(D) == 1) D <- diag(rep(1/n, n), nrow = n) if(length(D) == n) D <- diag(D, nrow = n) if(length(D) == n2) D <- as.matrix(D) # transformation en -.5*d^2 if(dissim == F) { ma <- max(X) diag(X) <- ma X <- ma - X } X <- (-0.5 * X * X) if(impres) { if(dissim == F) { cat(" - Positionement Multidimensionnel d'une matrice de similarite -\n" ) cat(" ---------------------------------------------------------------\n" ) cat("___________________________________________________________________________\n" ) } if(dissim == T) { cat(" - Positionement Multidimensionnel d'une matrice de dissimilarite -\n" ) cat(" ------------------------------------------------------------------\n" ) cat("______________________________________________________________________________\n" ) } } # initialisation valp <- NULL # vecteur des valeurs propres C <- NULL # matrice des composantes principales # calcul du W un <- matrix(1, nrow = n, ncol = 1) P <- diag(n) - un %*% t(un) %*% D W <- P %*% X %*% t(P) # decomposition de Choleski de D Dhalf <- chol(round(D, 6), pivot = T) Dhalf <- Dhalf[, order(attr(Dhalf, "pivot"))] Y <- Dhalf %*% W %*% Dhalf # diagonalisation eg <- eigen(Y, symmetric = T) valp <- eg$values vpos <- sum(valp > 1e-07) vnul <- sum(abs(valp) <= 1e-07) vneg <- sum(valp < (-1e-07)) inertiaX <- sum(valp[1:vpos]) # impression de l'histogramme des valeurs propres if(impres) { cat(paste("Inertie totale =", format(inertiaX), "\n")) cat("histogramme des valeurs propres (o ou n) ?\n") plth <- scan(quiet=T,"", character(), 1) if(length(plth) == 0) plth <- "n" if(plth == "o" || plth == "O") { par(mfrow = c(1, 1)) barplot(valp, space = 2, names = paste("v. p.", 1:n)) title("valeurs propres") } # choix du nombre d'axes si non donne en argument cat("__________________________________________________________________\n" ) valtab <- matrix(0, n, 3) dimnames(valtab) <- list(format(1:n), c("val.pro.", "% inert.", "% cumul.")) valtab[, 1] <- round(valp, digits = 4) valtab[, 2] <- round(valp/inertiaX * 100, digits = 2) for(i in 1:n) valtab[i, 3] <- sum(valtab[1:i, 2]) print(valtab) if(k == 0) { repeat { cat("Combien d'axes voulez-vous ? (<=", n, ")\n") k <- scan(quiet=T,"", n = 1) if(k != 0) break } } } # calcul des composantes principales Dhinv <- invgene(Dhalf) Cent <- Dhinv %*% eg$vectors[, 1:vpos] C <- Dhinv %*% eg$vectors[, 1:k] for(j in 1:k) C[, j] <- C[, j] * sqrt(valp[j]) for(j in 1:vpos) Cent[, j] <- Cent[, j] * sqrt(valp[j]) if(k == 1) C <- as.matrix(C) dimnames(C) <- list(dimnames(X)[[1]], paste("c", format(1:k), sep = "") ) # calcul et impression des aides a l'interpretation if(aideus == 1) { CTRus <- matrix(NA, nrow = n, ncol = k) COSus <- CTRus for(i in 1:n) for(j in 1:k) CTRus[i, j] <- round((D[i, i] * C[i, j]^2)/valp[ j] * 10000, digits = 0) dimnames(CTRus) <- list(dimnames(X)[[1]], paste("CTR", format(1: k), sep = "")) for(i in 1:n) { ss <- 0 for(j in 1:vpos) ss <- ss + Cent[i, j]^2 for(j in 1:k) COSus[i, j] <- round(Cent[i, j]^2/ss * 10000, digits = 0) } dimnames(COSus) <- list(dimnames(X)[[1]], paste("COS", format(1: k), sep = "")) } if((impres == T) | aideus == 1) { cat("_______________________________________________________________\n" ) cat("aides a l'interpretation pour les u.s. (o/n) ?\n") plta <- scan(quiet=T,"", character(), 1) if(length(plta) == 0) plta <- "n" if(plta == "o" || plth == "O") { cat("Contributions absolues des", n, "u.s. pour les", k, "premieres composantes\n") print(CTRus) cat("Contributions relative des", n, "u.s. pour les", k, "premieres composantes\n") print(COSus) } } # trace des eventuels graphiques if(graph == T) { cat("_______________________________________________________________\n" ) repeat { cat("graphique pour les u.s. (o/n) ?\n") pltc <- scan(quiet=T,"", character(), 1) if((length(pltc) == 0) | (pltc == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltch <- scan(quiet=T,"", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltcv <- scan(quiet=T,"", numeric(), 1) par(mfrow = c(1, 1), pty = "s") axespar <- c(pltch, pltcv) plot(C[, axespar], xlab = paste("c", axespar[1], " ", round(valp[axespar[1]], digits = 4), "(", round(valp[axespar[1]]/inertiaX * 100, digits = 2), "%)"), ylab = paste("c", axespar[2], " ", round(valp[axespar[2]], digits = 4), "(", round(valp[axespar[2]]/inertiaX * 100, digits = 2), "%)"), type = "n") abline(h = 0) abline(v = 0) text(C[, axespar], dimnames(C)[[1]]) } } # fin des graphiques } return(nomfichX, Xini, X, D, W, inertiaX, k, valp, C, CTRus, COSus) } ########################################## star.graph<-function(XY, Gr, D = 1, taille = 0.1) { # # Graphique en etoile # # Entrees # # XY matrice des coordonnees des points # Gr matrice des modalites d'appartenance des points # D metrique des variables, si D=1 alors 1/nIdn, si D est un vecteur alors D=diag(vecteur), sinon D=matrice # taille dimension en 'inches' des noms des classes sur les graphiques # # Sorties # Aucune # # Fonctions exterieures appelees # Aucune. # # R.Sabatier # MAJ 03/04/97 # nomfichXY <- deparse(substitute(XY)) XY <- as.matrix(XY) n <- nrow(XY) if(length(dimnames(XY)[[1]]) == 0) dimnames(XY)[[1]] <- paste("i", 1:n, sep = "") Gr <- as.matrix(Gr) q <- ncol(Gr) if(length(dimnames(Gr)[[2]]) == 0) dimnames(Gr)[[2]] <- paste("Var", 1:q, sep = "") n2 <- n * n k <- ncol(XY) if(k < 2) return() if(n != nrow(Gr)) return() # calcul de la metrique D if(length(D) == 1) D <- diag(rep(1/n, n), nrow = n) if(length(D) == n) D <- diag(D, nrow = n) if(length(D) == n2) D <- as.matrix(D) # verification du graphisme if(!exists(".Device", frame = 0)) { cat("Initialisez le graphique !!!\n") return() } # trace des graphiques { cat("_______________________________________________________________\n" ) repeat { cat("graphique pour les u.s. (o/n) ?\n") pltc <- scan(quiet=T,"", character(), 1) if(length(pltc) == 0)pltc<-"n" if((pltc == "n")) break else { # choix des numeros axes et de la variable qualitative cat("axe horizontal (<=", k, ") ?\n") pltch <- scan(quiet=T,"", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltcv <- scan(quiet=T,"", numeric(), 1) cat("numero de la var. qual. (<=", q, ") ?\n") pltvqual <- scan(quiet=T,"", numeric(), 1) # calcul des moyennes par classe nbmodal <- max(Gr[, pltvqual]) Cs <- matrix(nrow = nbmodal, ncol = k, 0) if(nbmodal <= 1) return() poicla <- matrix(ncol = 1, nrow = nbmodal, 0) for(i in 1:n) poicla[Gr[i, pltvqual]] <- poicla[Gr[i, pltvqual]] + D[i, i] for(i in 1:n) { cl <- Gr[i, pltvqual] poi <- D[i, i] x <- XY[i, pltch] y <- XY[i, pltcv] Cs[cl, pltch] <- Cs[cl, pltch] + x * poi Cs[cl, pltcv] <- Cs[cl, pltcv] + y * poi } for(i in 1:nbmodal) { no <- poicla[i] Cs[i, pltch] <- Cs[i, pltch]/no Cs[i, pltcv] <- Cs[i, pltcv]/no } dimnames(Cs) <- list(c(1:nbmodal), dimnames(XY)[[ 2]]) par(mfrow = c(1, 1), pty = "s") axespar <- c(pltch, pltcv) Ctot <- rbind(XY, Cs) plot(Ctot[, axespar], xlab = paste("c", axespar[ 1]), ylab = paste("c", axespar[2]), main = dimnames(Gr)[[2]][pltvqual], type = "p", pch = 18) for(i in 1:n) { cl <- Gr[i, pltvqual] segments(Cs[cl, pltch], Cs[cl, pltcv], XY[i, pltch], XY[i, pltcv], col = 7) } abline(h = 0) abline(v = 0) points(Cs[, axespar[1]], Cs[, axespar[2]], pch = 16, mkh = (taille + 0.09), col = 0) points(Cs[, axespar[1]], Cs[, axespar[2]], pch = 1, mkh = (taille + 0.09), col = 1) text(Cs[, axespar[1]], Cs[, axespar[2]], dimnames(Cs)[[1]], csi = taille, col = 3) } } } } ############################################## table.val<-function(N, fl = 1, fc = 1, ptypar = "s", type = 1, taille = 0.6, ...) { # # Graphique plan d'un tableau de contingence # # Entrees # # N tableau de contingence # type forme du graphique : 1 on represente le tableau de contingence # 2 on represente les profils lignes # 3 on represente les profils colonnes # 4 on represente nij/ni.n.j - 1 # taille facteur de proportionalite pour les graphiques # # R.Sabatier revu par JFD .... # lecture des donnees nomfichN <- deparse(substitute(N)) N <- as.matrix(N) # if(length(c(which.inf(N), which.na(N)))) # stop("valeurs manquantes ou infinies dans N") I <- nrow(N) J <- ncol(N) n <- sum(N) if(type < 1 | type > 4) return() # verification du graphisme if(!exists(".Device", frame = 0)) { cat("Initialisez le graphique !!!\n") return() } # graphique au carre ou rectangle par(pty = ptypar) # fichier si ordonancement selon fc et fl apres une analyse factorielle if(length(fl) != 1) { grfl <- 1 fl <- as.matrix(fl) if(nrow(fl) != I) return() ql <- ncol(fl) cat("numero de la var. ordre ligne (<=", ql, ") ?\n") pltfl <- scan(quiet=T,"", numeric(), 1) } else grfl <- 0 if(length(fc) != 1) { grfc <- 1 fc <- as.matrix(fc) if(nrow(fc) != J) return() qc <- ncol(fc) cat("numero de la var. ordre colonne (<=", qc, ") ?\n") pltfc <- scan(quiet=T,"", numeric(), 1) } else grfc <- 0 N <- N/n xy <- matrix(ncol = 2, nrow = I * J, NA) Z <- matrix(nrow = I * J, ncol = 1, NA) for(i in 1:I) for(j in 1:J) { ij <- i + (j - 1) * I if(grfc == 0) xy[ij, 1] <- j else xy[ij, 1] <- fc[j, pltfc] if(grfl == 0) xy[ij, 2] <- i else xy[ij, 2] <- fl[i, pltfl] if(type == 1) Z[ij] <- N[i, j] if(type == 3) Z[ij] <- N[i, j]/sum(N[, j]) if(type == 2) Z[ij] <- N[i, j]/sum(N[i, ]) if(type == 4) Z[ij] <- N[i, j]/(sum(N[, j]) * sum(N[i, ])) - 1 } # calcul de la taille des cercles et carres minx <- min(xy[, 1]) - 1.77 * taille * sqrt(abs(max(Z))) maxx <- max(xy[, 1]) + 1.77 * taille * sqrt(abs(max(Z))) miny <- min(xy[, 2]) - 1.77 * taille * sqrt(abs(max(Z))) maxy <- max(xy[, 2]) + 1.77 * taille * sqrt(abs(max(Z))) ming <- min(minx, miny) maxg <- max(maxx, maxy) # trace if(type == 2) { par(mar = c(5.1, 10, 4.1, 2.1), mgp = c(0, 5, 0)) plot(xy, type = "p", pch = 4, xlim = c(ming, maxg), ylim = c( ming, maxg), xlab = "", ylab = "", mkh = 0.01, bty = "l", axes = F, las = 1, main = "Profils", col = 1) axis(2, at = 1:I, labels = dimnames(N)[[1]], col = 1, ...) } if(type == 3) { par(mar = c(8, 4.1, 4.1, 2.1), mgp = c(0, 2, 0)) plot(xy, type = "p", pch = 4, xlim = c(ming, maxg), ylim = c( ming, maxg), xlab = "", ylab = "", mkh = 0.01, bty = "l", xaxt = "n", axes = F, las = 2, main = "Profils", col = 1) axis(1, at = 1:J, labels = format(1:J), col = 1, ...) } if((type == 1) | (type == 4)) { if(type == 1) titre <- "Table de Contingence" else titre <- "Tableau Pij / Pi.P.j - 1" plot(xy, type = "p", pch = 4, xlim = c(ming, maxg), ylim = c( ming, maxg), xlab = "", xaxt = "n", yaxt = "n", ylab = "", mkh = 0.01, bty = "l", main = titre, col = 1) axis(1, at = 1:J, labels = format(1:J), col = 1, ...) axis(2, at = 1:I, labels = dimnames(N)[[1]], col = 1, ...) } xyp <- matrix(nrow = I * J, ncol = 2, NA) xyn <- matrix(nrow = I * J, ncol = 2, NA) rp <- matrix(nrow = I * J, ncol = 1, NA) rn <- matrix(nrow = I * J, ncol = 1, NA) iip <- 1 iin <- 1 IJ <- I * J for(i in 1:IJ) { if(Z[i] > 0) { rp[iip] <- taille * sqrt(Z[i]) xyp[iip, ] <- xy[i, ] iip <- iip + 1 } if(Z[i] < 0) { rn[iin] <- 1.77 * taille * sqrt( - Z[i]) xyn[iin, ] <- xy[i, ] iin <- iin + 1 } } if(any(Z >= 0)) symbols(xyp[, 1], xyp[, 2], circles = rp, inches = F, add = T, xpd = T, smo = 0, col = 2) if(any(Z < 0)) symbols(xyn[, 1], xyn[, 2], squares = rn, inches = F, add = T, xpd = T) } ############################################# xy.val<-function(XY, taille = 1,...) { # # Graphique plan et representation d'une variable en z # # Entrees # # XY matrice des coordonnees des points # # R.Sabatier revu par JFD!!! # lecture des donnees nomfichXY <- deparse(substitute(XY)) if(!is.matrix(XY)) stop("XY n'est pas une matrice") XY <- as.matrix(XY) if(length(c(which.inf(XY), which.na(XY)))) stop("valeurs manquantes ou infinies dans XY") n <- nrow(XY) k <- ncol(XY) if(k < 3) return() # verification du graphisme if(!exists(".Device", frame = 0)) { cat("Initialisez le graphique !!!\n") return() } # trace des graphiques { cat("_______________________________________________________________\n" ) repeat { cat("graphique pour les u.s. (o/n) ?\n") pltc <- scan(quiet=T,"", character(), 1) if((length(pltc) == 0) | (pltc == "n")) break else { # choix des numeros axes et de la variable en z cat("axe horizontal (<=", k, ") ?\n") pltch <- scan(quiet=T,"", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltcv <- scan(quiet=T,"", numeric(), 1) cat("numero de la var. en Z (<=", k, ") ?\n") pltcz <- scan(quiet=T,"", numeric(), 1) # calcul de la taille des cercles et carres minxXY <- min(XY[, pltch]) maxxXY <- max(XY[, pltch]) minyXY <- min(XY[, pltcv]) maxyXY <- max(XY[, pltcv]) minzXY <- min(XY[, pltcz]) maxzXY <- max(XY[, pltcz]) minxp <- minxXY - taille * abs(max(XY[, pltcz]) ) maxxp <- maxxXY + taille * abs(max(XY[, pltcz]) ) minyp <- minyXY - taille * abs(max(XY[, pltcz]) ) maxyp <- maxyXY + taille * abs(max(XY[, pltcz]) ) # trace axespar <- c(pltch, pltcv) par(mfrow=c(1,1),pty="s") plot(XY[, axespar], xlab = paste("c", axespar[1 ]), ylab = paste("c", axespar[2]), type = "p", pch = 4, xlim = c(minxp, maxxp), ylim = c( minyp, maxyp), type = "n",col=1) abline(h = 0,col=1) abline(v = 0,col=1) r <- XY[, pltcz] XYp <- matrix(nrow = n, ncol = 2, NA) XYn <- matrix(nrow = n, ncol = 2, NA) rp <- matrix(nrow = n, ncol = 1, NA) rn <- matrix(nrow = n, ncol = 1, NA) iip <- 1 iin <- 1 for(i in 1:n) { if(r[i] > 0) { rp[iip] <- sqrt(r[i]) XYp[iip, 1] <- XY[i, pltch] XYp[iip, 2] <- XY[i, pltcv] iip <- iip + 1 } if(r[i] < 0) { rn[iin] <- sqrt( - r[i]) XYn[iin, 1] <- XY[i, pltch] XYn[iin, 2] <- XY[i, pltcv] iin <- iin + 1 } } if(any(XY[,pltcz]>=0)){ symbols(XYp[, 1], XYp[, 2], circles = taille * rp, inches = F, add = T, xpd = T, smo = 0, col = 2) text(XY[XY[, pltcz]>=0, axespar], dimnames(XY)[[1]][XY[, pltcz]>=0], col = 2,...) } if(any(XY[,pltcz]<0)){ symbols(XYn[, 1], XYn[, 2], squares = 1.77*taille * rn, inches = F, add = T, xpd = T, col=7 ) text(XY[XY[, pltcz]<0, axespar], dimnames(XY)[[1]][XY[, pltcz]<0], col = 7,...) } } } } } ################################################ Codisjc<-function(X) { # codage disjonctif complet de la matrice X # # Entree # X matrice (numerique) des modalites des variables qualitatives # Sorties # nbmodal vecteur des nombres de modalites de chaque variable # nbmodaltot nombre total de modalites # U matrice du codage # B matrice de Burt associee a U # X<-as.matrix(X) n<-nrow(X) p<-ncol(X) U<-NULL nbmodal<-NULL nommod<-NULL ylab<-NULL nbmodaltot<-as.numeric(0) labels<-rep(list(NULL),p) for(j in 1:p) { nbmodal[j]<-max(X[,j]) nommod<-paste(dimnames(X)[[2]][j],(1:nbmodal[j]),sep="") nbmodaltot<-nbmodaltot+nbmodal[j] V<-matrix(0,n,nbmodal[j]) for(i in 1:n) V[i,X[i,j]]<-1 U<-cbind(U,V) ylab<-c(ylab,nommod) labels[[j]]<-nommod } B<-t(U)%*%U nbmod<-nbmodal dimnames(U)<-list(dimnames(X)[[1]],ylab) dimnames(B)<-list(ylab,ylab) return(list(X=X,nbmod=nbmod,U=U,B=B,labels=labels)) } #################################################### discri<-function(X,g,Y,ResponseName="Y",imin=1,imax=nrow(Y),D=1, graph=T,ti=1,tj=2,cexpar=1,askpar=T) { # entrees: # X matrice dont les lignes sont les coordonnees de points(matrix of components) # g vecteur des modalites des groupes ou matrice du codage disjonctif complet # Y matrice des vecteurs lignes des coordonnees des points a affecter; ncol(Y)=ncol(X) # si Y manque, alors Y=X. # imin et imax indices des individus de Y a affecter. # D vecteur des poids des individus. # objectif: a quel groupe (a quelle modalite de g) affecter les lignes de Y ? # sortie : liste de composantes: # La matrice des centres de gravite des groupes, A. # La metrique de Mahalanobis de X, V. # Le vecteur des numero des groupes d'affectation des individus de Y, # groupe. # Si Y=X, (Y manquant) alors: # Le tableau croisant les groupes reels avec les groupes affectes, tableau. #-------------------------------------------------------------------- if(missing(Y)){ Y<-as.matrix(X) indicmiss=T } if(length(dim(g))>1) { if(dim(g)[2]==1) { indicator=rep(0,dim(g)[1]) indicator=g[,1] } if(sum(g)==dim(X)[1]) { indicator=rep(0,nrow(X)) for(i in 1:nrow(X)) for(j in 1:ncol(g)) if(g[i,j]==1)indicator[i]=j } g=indicator } # A matrice des coordonnees des centres de gravites # A est max(g) x ncol(X) ng<-max(g) codg<-codisj(g) G<-Mahalanobis(codg,D) if(length(D)==1) A<-G%*%t(codg)%*%X/nrow(X) else { XC<-D*X A<-G%*%t(codg)%*%XC } V<-Mahalanobis(X,D) #browser() groupe<-rep(0,nrow(Y)) for(j in imin:imax) { x<-Y[j,] U<-NULL for(i in 1:ng) U<-rbind(U,x) U<-(U-A) UU<-U%*%V aa<-NULL for(i in 1:ng) aa[i]<-sum(UU[i,]*U[i,]) distmin<-min(aa) i<-0 repeat {i<-i+1 if(aa[i]==distmin) break } groupe[j]<-i } tableau<-matrix(0,ng,ng) if(indicmiss) { for(i in 1:nrow(X)) tableau[groupe[i],g[i]]<-tableau[groupe[i],g[i]]+1 dimnames(tableau)<-list(paste("predicted",format(1:ng)),paste("real",format(1:ng))) if(graph) { cat("click on the plot to locate the legend\n") par(mfrow=c(1,1),pty="m",ask=askpar) if(length(ResponseName)==1) ResponseLevel=paste(ResponseName,1:max(g)) else {if(length(ResponseName)==max(g)) ResponseLevel=ResponseName else ResponseLevel=paste("Y",1:max(g)) } notgood=g!=groupe plot(X[,c(ti,tj)],type="n",xlab=dimnames(X)[[2]][ti],ylab=dimnames(X)[[2]][tj]) text(X[,c(ti,tj)],dimnames(X)[[1]],cex=cexpar,col=g+1) abline(h=0,v=0) points(A,pch=8,cex=cexpar*2,col=(min(g):max(g))+1) points(X[notgood,c(ti,tj)],pch=1,cex=cexpar*2.5) legend(locator(1),ResponseLevel,cex=cexpar,fill=(min(g):max(g))+1) } return(list(A=A,V=V,groupe=groupe,tableau=tableau)) } else return(list(A=A,V=V,groupe=groupe)) } #################################################### codisj<-function(x) { x<-as.vector(x) nbmodal<-max(x) X<-matrix(0,length(x),nbmodal) for(i in 1:length(x)) X[i,x[i]]<-1 return(X) } #################################################### MVcut<-function(X,qual,nbmod=3,equibreaks=F,breaks,labelsquant,labelsqual,graph=F,cexpar=1,labpar=c(5,7,7)) { # "M"ulti "V"ariables "cut"ing of quantitative variables # codage disjonctif complet de variables quantitatives et qualitatives placées en dernier #ENTREES # X matrice de variables quantitatives + d'éventuelles variables qualitatives (booleennes ou entières) #qual, entier numero de la première colonne qualitative (les autres suivent) #nbmod vecteur entiers, nombre de modalités de chaque variable quantitative #equibreaks vecteur booleen, si T séparateurs internes equidistants si F aux nbmod quantiles #breaks liste, liste des vecteurs de separateurs interieurs a chaque variable quantitative #labelsquant liste des noms des modalités des variables quantitatives #labelsqual liste des noms des modalités des variables qualitatives # labpar defaults to c(5,7,7), that is c(x,y,len) where x and y give the (approximate) number of tickmarks on the x and y axes # graph, T or F, if T, one can choose graphically the nb of levels and the location of break points when breaks is missing #SORTIES #X dataframe, tableau d'entiers des modalites #U matrice du codage disjonctif complet #nbmod vecteur du nb de modalités #breaks liste des vecteurs des séparateurs intérieurs + le max de la variable quantitavive #labels liste des noms des modalités des variables #EXEMPLE # # stateafcm<-MVcut(cbind(state.x77,Region=codes(state.region)),qual=ncol(state.x77)+1,labelsqual=sort(levels(state.region))) # ou bien # stateafcm<-MVcut(cbind(state.x77,Region=as.numeric(state.region)),qual=ncol(state.x77)+1,labelsqual=levels(state.region)) # afc(stateafcm,AFCM=T) # pour l'afcm du tableau des var quantitatives + la var qualitative des 4 modalités région # # par(pty="m",lab=labpar) X<-as.matrix(X) n<-nrow(X) pX<-ncol(X) if(!missing(qual)){ if(qual>1){ Xquant<-X[,1:(qual-1),drop=F] dimnames(Xquant)[[2]]=dimnames(X)[[2]][1:(qual-1)] } else { Xquant<-NULL U<-NULL } Xqual<-X[,qual:pX,drop=F] } else { Xquant<-X Xqual<-NULL } if(is.matrix(Xquant)) {# p<-ncol(Xquant) if(length(nbmod)==1) nbmod<-rep(nbmod,p) if(length(equibreaks)==1)equibreaks<-rep(equibreaks,p) if(missing(breaks)) { breaks<-rep(list(NULL),p) for(i in 1:p) { if(equibreaks[i]) {etend<-range(Xquant[,i]) hh<-diff(etend)/nbmod[i] for(j in 1:(nbmod[i])) breaks[[i]][j]<-etend[1]+j*hh } else { breaks[[i]]<-quantile(Xquant[,i],probs=(1:nbmod[i])/nbmod[i]) } endbreaks=breaks[[i]][nbmod[i]] if(graph) { titre="breaks : " for(j in 1:(nbmod[i]-1)) titre=c(titre,as.character(round(breaks[[i]][j],2))) titi=cut(Xquant[,i],breaks=c(min(Xquant[,i])-1,breaks[[i]]),labels=F) plot(Xquant[,i],ylab=dimnames(X)[[2]][i],type="n") text(Xquant[,i],dimnames(Xquant)[[1]],cex=cexpar,col=titi+1) abline(h=breaks[[i]][1:(nbmod[i]-1)],lty=6,lwd=2) title(main=titre) repeat{ cat(paste(dimnames(X)[[2]][i]," :\n")) cat(paste("nunber of levels : ",nbmod[i],"\n")) cat("Change or not that number ? (y/n)") repo1<-scan(quiet=T,what=character(),n=1) if(length(repo1)==0) repo1="n" if((repo1=="y")|(repo1=="Y")) {#y1 cat("nb of levels for ",dimnames(X)[[2]][i], " : ") nbmod[i]=scan(quiet=T,what=numeric(),n=1) cat("\n") }#y1 end if(repo1=="n") {cat("Try new",nbmod[i]-1 ,"break points for ",dimnames(X)[[2]][i],"(y,Y/n,N)?") repo<-scan(quiet=T,what=character(),n=1) } else { if(nbmod[i]==2) cat("Try new",nbmod[i]-1 ,"break point for ",dimnames(X)[[2]][i]) else cat("Try new",nbmod[i]-1 ,"break points for ",dimnames(X)[[2]][i],"\n") repo="y" } if(length(repo)==0)repo="n" if((repo=="y")|(repo=="Y")) {#y # repeat{ if(nbmod[i]>2)cat("the breaks :\n") toto<-scan(quiet=T,"", numeric(),n=nbmod[i]) #cat("OK (y,Y/n,N) for ",toto," ") #repon<-scan(quiet=T,what=character(),n=1) #if(length(repon)==0)repon<-"y" #if((repon=="y")|(repon=="Y")) break #if(length(toto)==nbmod[i])break # } breaks[[i]]<-c(toto,endbreaks) titre="breaks : " for(j in 1:(nbmod[i]-1)) titre=c(titre,as.character(round(breaks[[i]][j],2))) titi=cut(Xquant[,i],breaks=c(min(Xquant[,i])-1,breaks[[i]]),labels=F) plot(Xquant[,i],ylab=dimnames(X)[[2]][i],type="n") text(Xquant[,i],dimnames(Xquant)[[1]],cex=cexpar,col=titi+1) abline(h=breaks[[i]][1:(nbmod[i]-1)],lty=4,lwd=2) title(titre) #browser() }#end y if(repo=="n") { cat("------------------------------\n") break } }#end repeat }#end graph }# end for }#end missing breaks else {nbmod<-rep(0,p) for(i in 1:p) {breaks[[i]]<-c(breaks[[i]],max(Xquant[,i])) nbmod[i]<-length(breaks[[i]]) } } #browser() if(missing(labelsquant)){ labelsquant<-rep(list(NULL),p) for(i in 1:p) labelsquant[[i]]<-paste(dimnames(Xquant)[[2]][i],format(1:nbmod[i]),sep="") } mat<-Xquant for(i in 1:p) mat[,i]<-cut(Xquant[,i],breaks=c(min(Xquant[,i])-1,breaks[[i]]),labels=F) codage<-Codisjc(mat) X<-mat U<-codage$U }# if(is.matrix(Xqual)) {## cat("variables Entières ou Codage disjonctif complet?(E ou C)\n") pltc <- scan(quiet=T,"", character(), 1) if((pltc=="E")|(pltc=="e")) {codage1<-Codisjc(Xqual) V<-codage1$U if(missing(labelsqual)) labelsqual<-codage1$labels else dimnames(V)[[2]]<-labelsqual U<-cbind(U,V) nbmodqual<-codage1$nbmod } else { pqual<-sum(Xqual[1,]) if(pqual==1)cat("rentrez le nombre de modalités de cette variable\n") else cat("Rentrez le nombre de modalités pour chacune des",pqual,"variables\n") nbmodqual<-scan(quiet=T,"",numeric(),pqual) U<-cbind(U,Xqual) if(missing(labelsqual)) labelsqual<-dimnames(Xqual)[[2]] } }## #B<-t(U)%*%U if((is.matrix(Xquant))&(is.matrix(Xqual))) { if((pltc=="E")|(pltc=="e")) X<-cbind(mat,Xqual) labels<-c(labelsquant,labelsqual) nbmod<-c(nbmod,nbmodqual) } if((!is.matrix(Xquant))&(is.matrix(Xqual))) { labels<-labelsqual nbmod<-nbmodqual } if((is.matrix(Xquant))&(!is.matrix(Xqual))) labels<-labelsquant X<-as.data.frame(X) if(is.matrix(Xquant)) { for(i in 1:dim(Xquant)[2]) for(j in 1:nbmod[i]) { if(sum(U[,j])==0) { cat("the variable ",dimnames(Xquant)[[2]][i]," is of ",j,"th empty level!\n") return() } } return(list(X=X,U=U,nbmod=nbmod,labels=labels,breaks=breaks)) } else return(list(U=U,nbmod=nbmod,labels=labels)) } #################################################### myps<-function(file,largeur=6,hauteur=6,horizpar=F,devicepar=windows) { # création d'un fichier postscript d'une image R # #file chaine de caracteres, nom du fichier postscript #devicepar chaine de caractères, nom du device de l'image source, windows ou X11 # # Pour visualiser et/ou imprimer le fichier postscript utiliser ghostview # création de J.F. Durand, le 17/03/2001 # dev.copy(device=devicepar) dev.print(device=postscript,file=paste(file,".ps",sep=""),width=largeur,height=hauteur,horizontal=horizpar) dev.off(dev.prev()) } #################################################### mygraphics<-function(file,width=6,height=6,horizontal=F,devicepar=X11) { # création d'un fichier postscript d'une image R # #file chaine de caracteres, nom du fichier postscript, ou pdf ou jpeg #devicepar chaine de caractères, nom du device de l'image source, windows ou X11 #horizontal, logical, uniquement en postscript, T donne paysage, F portrait # # Pour visualiser et/ou imprimer le fichier postscript utiliser ghostview # création de J.F. Durand, le 30/09/2004 # cat("Dimensions : width =",width,", height =",height,"\n") reponse<-menu(c(" postscript (dimensions in inches, 1 inch ~ 2.54 cm)", " pdf (dimensions in inches, 1 inch ~ 2.54 cm)", " jpeg (dimensions in pixels)"),title=paste("Choose the format of the file '",file,"' ( 0 to exit ):",sep="")) if(reponse==0)invisible(return()) if(reponse==1) { dev.copy(device=devicepar) dev.print(device=postscript,file=paste(file,".ps",sep=""),width=width,height=height, horizontal=horizontal,bg=par("bg")) } if(reponse==2) { dev.copy(device=devicepar) dev.print(device=pdf,file=paste(file,".pdf",sep=""),width=width,height=height,bg=par("bg")) } if(reponse==3) { dev.copy(device=devicepar) dev.print(device=jpeg,file=paste(file,".jpeg",sep=""),width=width,height=height,bg=par("bg")) } dev.off(dev.cur()) } #################################################### Profils<-function(TC,cexpar=1,titlepar=T,bgpar="lavender") { # édition des diagrammes en batons des profils et des marges d'un tableau de contingence # TC, matrice, tableau de contingence # cexpar, scalaire, paramètre cex pour l'édition des graphiques # titlepar, booléen, si T le titre du grahique est édité. # bgpar couleur du background graphique, utiliser aussi "cornsilk" ou "lightgray" utiliser aussi colours()[624] # # Auteur : Jean-Francois Durand 10/10/2001 # TC<-as.matrix(TC) I<-nrow(TC) J<-ncol(TC) N<-sum(TC) P<-TC/N margeC<-apply(P,2,sum) margeL<-apply(P,1,sum) ProfilsL<-P/margeL ProfilsC<-t(t(P)/margeC) par(ask=T) par(bg=bgpar) palette(rainbow(12)) cat("Diagrammes des Profils Lignes + diagramme marge Colonnes\n") cat("Nombre de graphes par ligne ? (<=",I+1,")\n") colo<-scan(quiet=T,"",numeric(),1) cat("Nombre de lignes ?\n") lign<-scan(quiet=T,"",numeric(),1) par(oma=c(0,0,4,0),mfrow=c(lign,colo),pty="s") for(i in 1:I) barplot(ProfilsL[i,],ylim=c(0,max(rbind(ProfilsL,margeC))),main=paste(dimnames(TC)[[1]][i]," (",round(margeL[i],3),")",sep=""),cex=cexpar) barplot(margeC,ylim=c(0,max(rbind(ProfilsL,margeC))),main="Profil Moyen",cex=cexpar) if(titlepar) mtext(side=3,line=0,cex=cexpar+0.5,outer=TRUE,paste("Profils Lignes ","+ Marge Colonne",sep="")) cat("Diagrammes des Profils Colonnes + diagramme marge Ligne\n") cat("Nombre de graphes par ligne ? (<=",J+1,")\n") colo<-scan(quiet=T,"",numeric(),1) cat("Nombre de lignes ?\n") lign<-scan(quiet=T,"",numeric(),1) par(oma=c(0,0,4,0),mfrow=c(lign,colo),pty="s") for(i in 1:J) barplot(ProfilsC[,i],ylim=c(0,max(cbind(ProfilsC,margeL))),main=paste(dimnames(TC)[[2]][i]," (",round(margeC[i],3),")",sep=""),cex=cexpar) barplot(margeL,ylim=c(0,max(cbind(ProfilsC,margeL))),main="Profil Moyen",cex=cexpar) if(titlepar) mtext(side=3,line=0,cex=cexpar+0.5,outer=TRUE,paste("Profils Colonnes ","+ Marge Ligne",sep="")) return(list(P=P,margeL=margeL,margeC=margeC,ProfilsL=ProfilsL,ProfilsC=ProfilsC)) } #################################################### PCA<-function(X, centrer = T, cor = T,k=3, impres = T, graph = T, Xl = 1, Xc = 1, aideus = 1, aideva = 1,colpar=0,bgpar="whitesmoke",askpar=F,cexpar=0.8) { # # usual Principal Component Analysis # # Entrees # # X matrice des variables # D metrique des variables, si D=1 alors 1/nIdn, si D est un vecteur alors D=diag(vecteur), sinon D=matrice # centrer=F on ne centre pas X, sinon centrer=T # cor=F on ne reduit pas, centrage et reduction si cor=T # k nombre de composantes retenues, par defaut 3 # impres impression des resultats si T, si F rien # graph trace des graphiques si T, si F pas de trace # Xl matrice eventuelle des u.s. supplememtaires # Xc matrice eventuelle des variables supplememtaires # aideus=1 aides a l'interpretation pour les u.s., CTR et COS, pour les k premiers axes, sinon aideus=0 # aideva=1 aides a l'interpretation pour les variables, CTR et COS, pour les k premiers axes, sinon aideva=0 # bgpar chaine de caractères, couleur du background, par exple, "bisque1", "whitesmoke"...utiliser colors() pour choisir # Sorties # # nomfichX nom du fichier X # Xini tableau initial # k rang de la DSD retenue # valp vecteur de toutes les valeurs propres # A matrice des k premiers axes principaux de la DSD # C matrice des k premieres composantes principales de la DSD # CRTus et COSus contributions absolues et relatives des u.s. pour les k premieres composantes # CTRva et COSva contributions absolues et relatives des variables pour les k premiers axes # As matrice des k premiers axes pour les eventuelles variables supplementaires # Cs matrice des k premieres composantes pour les eventuelles u.s. supplementaires # # Fonctions exterieures appelees # # Dcentred # # R.Sabatier revu par JF Durand # MAJ 03/04/97 20/01/2002 # lecture des donnees nomfichX <- deparse(substitute(X)) X <- as.matrix(X) # if(length(c(which.inf(X), which.na(X)))) # stop("valeurs manquantes ou infinies dans X") if(impres){ oldpar<-par(no.readonly = TRUE) par(ask=F,bg=bgpar,mfrow=c(1,1),bty="n",xaxt="n",yaxt="n") plot(0,0,type="n",bg=par("bg"),xlab="",ylab="",xlim=c(-4,4),ylim=c(-4,4)) text(0,2,"Principal Component",cex=3,col="red") text(0,0,"Analysis",cex=3,col="red") text(0,-2,"Jean-François Durand",cex=1.3,col="black") text(0,-3,"http://www.jf-durand-pls.com",cex=1,col="black") par(oldpar) par(ask=askpar,bg=bgpar,xaxt="s",yaxt="s") } D=1 Q=1 n <- nrow(X) p <- ncol(X) if(is.null(dimnames(X))) dimnames(X) <- list(format(1:n), paste("v", 1:p, sep = "")) if(length(dimnames(X)[[1]]) == 0) dimnames(X)[[1]] <- format(1:n) if(length(dimnames(X)[[2]]) == 0) dimnames(X)[[2]] <- paste("v", 1:p, sep = "") Xini <- as.matrix(X) if(p < 2) return() if(impres == F) graph <- F if(impres == T) { if(!exists(".Device", frame = 0)) { cat("Initialisez le graphique !!!\n") return() } } if(centrer == F) cor <- F if(cor == T) centrer <- T p2 <- p * p n2 <- n * n As <- NULL Cs <- NULL CTRus <- NULL COSus <- NULL CTRva <- NULL COSva <- NULL # calcul de la metrique Q if(length(Q) == 1) Q <- diag(p) if(length(Q) == p) Q <- diag(Q, nrow = p) if(length(Q) == p2) Q <- as.matrix(Q) # calcul de la metrique D if(length(D) == 1) D <- diag(rep(1/n, n), nrow = n) if(length(D) == n) D <- diag(D, nrow = n) if(length(D) == n2) D <- as.matrix(D) # verifications u.s. et variables supplementaires if(length(Xl) != 1) { isup <- 1 Xl <- as.matrix(Xl) if(ncol(Xl) != ncol(X)) return() if(is.null(dimnames(Xl))) dimnames(Xl) <- list(paste("is", 1:n, sep = ""), paste( "v", 1:p, sep = "")) if(length(dimnames(Xl)[[1]]) == 0) dimnames(Xl)[[1]] <- paste("is", 1:n, sep = "") if(length(dimnames(Xl)[[2]]) == 0) dimnames(Xl)[[2]] <- paste("v", 1:p, sep = "") } else isup <- 0 if(length(Xc) != 1) { vsup <- 1 Xc <- as.matrix(Xc) if(nrow(Xc) != nrow(X)) return() #browser() if(is.null(dimnames(Xc))) dimnames(Xc) <- list(paste("i", 1:n, sep = ""), paste( "vs", dim(Xc)[2], sep = "")) if(length(dimnames(Xc)[[1]]) == 0) dimnames(Xc)[[1]] <- paste("i", 1:n, sep = "") if(length(dimnames(Xc)[[2]]) == 0) dimnames(Xc)[[2]] <- paste("vs", 1:p, sep = "") } else vsup <- 0 if(centrer == T) { centrage <- Dcentred(X, D = D) X <- centrage$Xc meanX <- matrix(centrage$moy) varX <- matrix(centrage$var) dimnames(meanX) <- list(dimnames(X)[[2]], "moy") dimnames(varX) <- list(dimnames(X)[[2]], "var") if(isup == 1) Xl <- sweep(Xl, 2, meanX) if(vsup == 1) { centragevs <- Dcentred(Xc, D = D) Xc <- centragevs$Xc } } if(cor == T) { X <- centrage$Xcr if(isup == 1) Xl <- sweep(Xl, 2, sqrt(varX), FUN = "/") if(vsup == 1) Xc <- centragevs$Xcr } if(impres) { cat("__________________________________________________________________\n") cat("\n") cat(" - Principal Component Analysis -\n") cat(" -------------------------------\n") if(isup == 1) cat(" with supplementary observations\n") if(vsup == 1) cat(" with supplementary variables\n") cat("__________________________________________________________________\n" ) if(cor == T) { cat(" standardized data\n") cat(" mean of the variables\n") print(t(meanX)) cat(" variance of the variables\n") print(t(varX)) } if(centrer == T && cor == F) { cat(" centered data\n") cat(" men of the variables\n") print(t(meanX)) cat(" variance of the variables\n") print(t(varX)) } cat("__________________________________________________________________\n" ) } # initialisation valp <- NULL # vecteur des valeurs propres A <- NULL # matrice des facteurs principaux C <- NULL # matrice des composantes principales # decomposition de Choleski de Q Qhalf <- chol(round(Q, 6)) # Qhalf <- chol(round(Q, 6), pivot = T) # Qhalf <- Qhalf[, order(attr(Qhalf, "pivot"))] Y <- X %*% t(Qhalf) # diagonalisation diago <- t(Y) %*% D %*% Y eg <- eigen(diago, symmetric = T) valp <- eg$values if(all(eg$vectors[,1]<=0)) eg$vectors[,1]=-eg$vectors[,1] for(i in 1:ncol(diago)) if(valp[i] < 0) valp[i] <- 0 inertiaX <- sum(valp) valp2 <- sum(valp^2) # impression de l'histogramme des valeurs propres if(impres) { par(bg=bgpar) cat(paste("Total Inertia =", format(inertiaX), "\n" )) # choix du nombre d'axes si non donne en argument cat("__________________________________________________________________\n" ) valtab <- matrix(0, p, 3) dimnames(valtab) <- list(format(1:p), c("eigenvalues", " % inertia", " % cum. inertia")) valtab[, 1] <- round(valp, digits = 4) valtab[, 2] <- round(valp/inertiaX * 100, digits = 2) valtab[, 3] <- round(cumsum(valp/inertiaX * 100),2) # for(i in 1:p) { #valtab[i, 3] <- sum(valtab[1:i, 2]) # valtab[i, 4] <- sqrt(sum(valp[1:i]^2)/valp2) # calcul du RV #} print(valtab) cat("barplots of the eigenvalues (y or n) ?\n") plth <- scan(quiet=T,"", character(), 1) if(length(plth) == 0) plth <- "n" if(plth == "y" || plth == "Y") { par(mfrow = c(1, 2),pty="m") barplot(valp, space = 2, names = format(1:p),bg=par("bg"),cex.names =cexpar,cex.axis=cexpar,col=palette()[1:length(valp)+1]) title("eigenvalues",cex.main=cexpar+0.4) barplot(valtab[,3], space = 2, names = format(1:p),bg=par("bg"),cex.names =cexpar,cex.axis=cexpar,col=palette()[1:length(valp)+1]) abline(h=100,lty=3,col="red") title("% cumul. inertia",cex.main=cexpar+0.4) } repeat { cat("How many Components? (<=", p, ")") k <- scan(quiet=T,"", n = 1) if(k != 0) break } } # calcul des composantes principales Cent <- Y %*% eg$vectors C <- Cent[, 1:k] if(k == 1) C <- as.matrix(C) dimnames(C) <- list(dimnames(X)[[1]], paste("c", format(1:k), sep = "") ) # calcul des axes principaux Aent <- matrix(NA, nrow = ncol(X), ncol = p) if(length(Q) != 1) fac <- solve(Qhalf) %*% eg$vectors else fac <- eg$vectors for(i in 1:p) Aent[, i] <- (fac[, i]) * sqrt(valp[i]) A <- Aent[, 1:k] dimnames(A) <- list(dimnames(X)[[2]], paste("a", format(1:k), sep = "") ) # calcul des facteurs principaux Fa <- matrix(NA, nrow = ncol(X), ncol = k) Fa <- Q %*% A dimnames(Fa) <- list(dimnames(X)[[2]], paste("f", format(1:k), sep = "" )) # calcul des coordonnees u.s. et variables supplementaires if(isup == 1) { Cs <- Xl %*% Q %*% A for(i in 1:k) Cs[, i] <- (Cs[, i])/sqrt(valp[i]) dimnames(Cs) <- list(dimnames(Xl)[[1]], paste("c", format(1:k), sep = "")) } if(vsup == 1) { As <- t(Xc) %*% D %*% C for(i in 1:k) As[, i] <- (As[, i])/sqrt(valp[i]) dimnames(As) <- list(dimnames(Xc)[[2]], paste("a", format(1:k), sep = "")) } # calcul et impression des aides a l'interpretation if(aideus == 1) { CTRus <- matrix(NA, nrow = n, ncol = k) COSus <- CTRus for(i in 1:n) for(j in 1:k) CTRus[i, j] <- round((D[i, i] * C[i, j]^2)/valp[ j] * 10000, digits = 0) dimnames(CTRus) <- list(dimnames(X)[[1]], paste("CTA", format(1: k), sep = "")) for(i in 1:n) { ss <- 0 for(j in 1:p) ss <- ss + Cent[i, j]^2 for(j in 1:k) COSus[i, j] <- round(Cent[i, j]^2/ss * 10000, digits = 0) } dimnames(COSus) <- list(dimnames(X)[[1]], paste("COS", format(1: k), sep = "")) } if(aideva == 1) { CTRva <- matrix(NA, nrow = p, ncol = k) COSva <- CTRva for(i in 1:p) for(j in 1:k) CTRva[i, j] <- round((Q[i, i] * A[i, j]^2)/valp[ j] * 10000, digits = 0) dimnames(CTRva) <- list(dimnames(X)[[2]], paste("CTA", format(1: k), sep = "")) for(i in 1:p) { ss <- 0 for(j in 1:p) ss <- ss + Aent[i, j]^2 for(j in 1:k) COSva[i, j] <- round(Aent[i, j]^2/ss * 10000, digits = 0) } dimnames(COSva) <- list(dimnames(X)[[2]], paste("COS", format(1: k), sep = "")) } if(impres) { if(aideus == 1) { cat("_______________________________________________________________\n" ) cat("informations about the projected observations (y/n) ?\n") plta <- scan(quiet=T,"", character(), 1) if(length(plta) == 0) plta <- "n" if(plta == "y" || plta == "Y") { cat("Absolute Contributions of", n, " observations on the first ", k, "components\n" ) print(CTRus) cat("Relative Contributions of the", n, "observations on the first ", k, "components\n" ) print(COSus) } } if(aideva == 1) { cat("_______________________________________________________________\n" ) cat("informations about the projected variables (y/n) ?\n" ) plta <- scan(quiet=T,"", character(), 1) if(length(plta) == 0) plta <- "n" if(plta == "y" || plta == "Y") { cat("Absolute Contributions of the", p, "variables on the first", k, "axes\n") print(CTRva) cat("Relative Contributions relative of the ", p, "variables on the first", k, "axes\n") print(COSva) } } } # trace des eventuels graphiques if(graph == T) { cat("_______________________________________________________________\n" ) repeat { cat("scatterplots of the observations (y/n) ?\n") pltc <- scan(quiet=T,"", character(), 1) if(length(pltc) == 0) pltc<-"n" if((pltc == "N") | (pltc == "n")) break else { cat("horizontal axis (<=", k, ") ?\n") pltch <- scan(quiet=T,"", numeric(), 1) cat("vertical axis (<=", k, ") ?\n") pltcv <- scan(quiet=T,"", numeric(), 1) par(mfrow = c(1, 1), pty = "s") axespar <- c(pltch, pltcv) if(isup == 1) Ctot <- rbind(C, Cs) else Ctot <- C par(pty="m") plot(Ctot[, axespar], xlab = paste("t", axespar[ 1], " ", round(valp[axespar[1]], digits = 4), "(", round(valp[axespar[1]]/inertiaX * 100, digits = 2), "%)"), ylab = paste("t", axespar[ 2], " ", round(valp[axespar[2]], digits = 4), "(", round(valp[axespar[2]]/inertiaX * 100, digits = 2), "%)"), type = "n",bg=par("bg"),cex=cexpar) abline(h = 0) abline(v = 0) text(Ctot[, axespar], dimnames(Ctot)[[1]],col=colpar+1,cex=cexpar) cat("Map of the observations colored according to the levels of a categorical variable? (y/n)\n") plto <- scan(quiet=T,"", character(), 1) if( (length(plto)!=0)&&((plto=="Y")|(plto=="y"))) { repeat{ cat("Enter the name of the vector or that of the matrix containing the column-variable\n") nom <- scan(quiet=T,"", character(), 1) matricenom<-as.matrix(get(nom,pos=1)) if(nrow(matricenom)==nrow(Xini))break else cat("misfit on the number of rows\n") } if(is.vector(get(nom,pos=1))) {dimnames(matricenom)<-list(dimnames(Xini)[[1]],c(nom)) numer<-1 variab<-get(nom,pos=1) } if((is.matrix(get(nom,pos=1)))|(is.data.frame(get(nom,pos=1)))) { cat(paste(dimnames(matricenom)[[2]],"(",1:ncol(matricenom),"),",sep=""),"\n") cat("enter the column number of integer values\n") numer <- scan(quiet=T,"",numeric(), 1) variab<-matricenom[,numer] cat("name of the variable : ",dimnames(matricenom)[[2]][numer],"\n") } cat("click to locate the top left corner of the legend\n") for(i in min(variab):max(variab)) text(C[variab==i,axespar],dimnames(C)[[1]][variab==i],col=i+1,cex=cexpar) legend(locator(1),paste(dimnames(matricenom)[[2]][numer],min(variab):max(variab),sep=""), fill=(min(variab):max(variab))+1,cex=cexpar) } } } cat("_______________________________________________________________\n" ) repeat { if(cor == F) { cat("plots for the variables (y/n) ?\n") plta <- scan(quiet=T,"", character(), 1) if(length(plta) == 0) plta<-"n" if((plta == "N") | (plta == "n")) break else { cat("horizontal axis(<=", k, ") ?\n") pltah <- scan(quiet=T,"", numeric(), 1) cat("vertical axix (<=", k, ") ?\n") pltav <- scan(quiet=T,"", numeric(), 1) par(mfrow = c(1, 1), pty = "s") axespar <- c(pltah, pltav) if(vsup == 1) Atot <- rbind(A, As) else Atot <- A plot(rbind(Atot[, axespar],c(0,0)), xlab = paste("t", axespar[1], " ", round(valp[axespar[1]], digits = 4), "(", round(valp[axespar[1]]/ inertiaX * 100, digits = 2), "%)"), ylab = paste("t", axespar[2], " ", round(valp[ axespar[2]], digits = 4), "(", round(valp[ axespar[2]]/inertiaX * 100, digits = 2), "%)"), type = "n",bg=par("bg"),cex=cexpar) arrows(rep(0,p),rep(0,p),0.92*Atot[,axespar[1]],0.92*Atot[,axespar[2]],col=1,angle=15,length=0.15) abline(v=0) abline(h=0) text(Atot[, axespar], dimnames(Atot)[[1]],col=colpar+1,cex=cexpar) if(vsup==1){ #arrows(rep(0,nrow(As)),rep(0,nrow(As)),0.92*As[,axespar[1]],0.92*As[,axespar[2]],col=colpar+2,angle=15,length=0.15) text(As[, axespar,drop=F], dimnames(Xc)[[2]],col=colpar+2) } } } else { cat("circle of correlations (y/n) ?\n") plta <- scan(quiet=T,"", character(), 1) if(length(plta) == 0) plta<-"n" if((plta == "N") | (plta == "n")) break else { cat("horizontal axis (<=", k, ") ?\n") pltah <- scan(quiet=T,"", numeric(), 1) cat("vertical axix(<=", k, ") ?\n") pltav <- scan(quiet=T,"", numeric(), 1) par(mfrow = c(1, 1), pty = "s") axespar <- c(pltah, pltav) theta <- seq(0, 20, 0.05) x <- cos(theta) y <- sin(theta) if(vsup == 1) Atot <- rbind(A, As) else Atot <- A plot(x, y, type = "l", xlab = paste("t", axespar[1], " ", round(valp[axespar[1]], digits = 4), "(", round(valp[axespar[1]]/ inertiaX * 100, digits = 2), "%)"), ylab = paste("t", axespar[2], " ", round(valp[ axespar[2]], digits = 4), "(", round(valp[ axespar[2]]/inertiaX * 100, digits = 2), "%)"),cex=cexpar) abline(h = 0) abline(v = 0) arrows(rep(0,p),rep(0,p),0.92*Atot[, axespar[1]],0.92*Atot[, axespar[2]],col=1,angle=15,length=0.15) text(Atot[, axespar], dimnames(Atot)[[1]],col=colpar+1,cex=cexpar) if(vsup==1){ #arrows(rep(0,nrow(As)),rep(0,nrow(As)),0.92*As[,axespar[1]],0.92*As[,axespar[2]],col=colpar+2,angle=15,length=0.15) text(As[, axespar,drop=F], dimnames(Xc)[[2]],col=colpar+2,cex=cexpar) } } } } # fin des graphiques } invisible(return(list(nomfichX=nomfichX,Xini=Xini,k=k,valp=valp,A=A,Components=C,CTRus=CTRus, COSus=COSus,CTRva=CTRva,COSva=COSva,As=As,Cs=Cs))) } ########################################################## Hotelling<-function(X,H=ncol(X),FSthreshold=0.95,ellipse=T,ii=1,jj=2,nbpoints=100,cexpar=0.8,titlepar=T) { ###### #INPUTS #X is the matrix of components #H the selected optimal number of components (optimal PRESS-dimension),defaults to ncol(X) #FSthreshold is the Fisher-Snedecor confidence level, defaults to 0.95 #ellipse, boolean, T for displaying the FSthreshold-Hotelling ellipse on the t(ii),t(jj) scatterplot #nbpoints defaults to 100 used to draw the ellipse with 100 points. #ii and jj, integers (smaller than H) indicating what ii,jj pair of components to display, default to 1 and 2. # OUTPUTS #T2 is the n by H matrix of the Hotelling T^2 of the n observations computed using components 1, 1&2,...,1&2&...&H #threshold is the vector of length H giving the F-S thresholds when using 1,2,...,H components #outliers is a 0/1 vector of length n, 1 indicating that the observation is an H-outlier ###### X=as.matrix(X) n<-nrow(X) outliers<-rep(0,n) # outliers calculated on H components, if a value is 1 the observation is an H-outlier T2<-matrix(0,nrow=n,ncol=H) dimnames(T2)=list(dimnames(X)[[1]],seq(1,H)) threshold=rep(0,H) #Hotelling's T2 statistic based on H components #### for(i in 1:H){ S<-var(X[,1:i])# variance, division by n-1 Sinv<-solve(as.matrix(S)) XSX<-(n/(n-1))*X[,1:i]%*%Sinv%*%t(X[,1:i]) T2[,i]<-diag(XSX) threshold[i]<-(i*(n^2-1)/(n*(n-i)))*qf(FSthreshold,i,n-i,lower.tail=T) # Obtain upper threshold percentile for Hotelling's T2 based on selected components }#end for #### #comparison with F critical: H-outliers for(i in 1:n){ if((T2[i,H])>=threshold[H]) outliers[i]<-1 }# end for # Detecting and printing eventual H-outliers ######## if(sum(outliers)>0) { out=(1:nrow(X))[outliers==1] print(paste("Outliers with",H,"components")) print(dimnames(X)[[1]][out]) cat("\n") } ######## #displaying the ellipse on the t(ii), t(jj) component's scatterplot if(ellipse) { par(pty="m",mfrow=c(1,1)) variance=apply(X,2,var)# variance, division by n-1 quadratic=variance[c(ii,jj)]*threshold[2] angles=seq(0, 2*pi, length.out=nbpoints) Xell=sqrt(quadratic[1])*cos(angles) Yell=sqrt(quadratic[2])*sin(angles) plot(c(X[,ii],Xell),c(X[,jj],Yell),type="n",xlab=paste("t",ii),ylab=paste("t",jj),cex=cexpar) abline(h=0,v=0) if(titlepar) title(main="Hotelling Ellipse T2 95%") points(Xell,Yell,type="l") text(X[,ii],X[,jj],dimnames(X)[[1]],cex=cexpar) }#end if ellipse ####### return(list(Components=X,T2=T2,H=H,threshold=threshold,FSthreshold=FSthreshold,outliers=outliers)) } ##########################################################