########################################################### # Liste des fonctions sous R (version 0.963): # 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 affecte MVcut myps Profils ########################################################### #.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(R,valp,alpha,vect,WDC,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) 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(Yhat,Yres,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(V,U,mean,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(X,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(Xc,Xcr,moy,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("", 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("", 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("", 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("", character(), 1) if((length(pltc) == 0) | (pltc == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltch <- scan("", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltcv <- scan("", 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("", character(), 1) if((length(plta) == 0) | (plta == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltah <- scan("", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltav <- scan("", 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("", character(), 1) if((length(plta) == 0) | (plta == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltah <- scan("", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltav <- scan("", 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("", character(), 1) if((length(plta) == 0) | (plta == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltah <- scan("", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltav <- scan("", 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(nomfichX, Xini, X, QX, nomfichY, Yini, Y, QY, D, inertiaACPVI, k, valp, A, C, 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=2,bgpar="lavender") { # # 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, utiliser colours()[624] par exemple # 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 # MAJ 03/04/97 # 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") 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 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" )) cat("barplot des valeurs propres (o ou n) ?\n") plth <- scan("", character(), 1) if(length(plth) == 0) plth <- "n" if(plth == "o" || plth == "O") { par(mfrow = c(1, 1)) barplot(valp, space = 2, names = format(1:p)) title("valeurs propres") } # 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) 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) if(k == 0) { repeat { cat("Combien d'axes voulez-vous ? (<=", p, ")\n") k <- scan("", 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("", 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("", 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("", character(), 1) if(length(pltc) == 0) pltc<-"n" if((pltc == "N") | (pltc == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltch <- scan("", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltcv <- scan("", 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("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(Ctot[, axespar], dimnames(Ctot)[[1]],col=colpar+3) cat("Une carte des individus actifs colories selon les modalites d'une variable qualitative? (o/n)\n") plto <- scan("", 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("", 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("",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) legend(locator(1),paste(dimnames(matricenom)[[2]][numer],min(variab):max(variab),sep=""), fill=(min(variab):max(variab))+1) } } } cat("_______________________________________________________________\n" ) repeat { if(cor == F) { cat("graphique pour les variables (o/n) ?\n") plta <- scan("", character(), 1) if(length(plta) == 0) plta<-"n" if((plta == "N") | (plta == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltah <- scan("", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltav <- scan("", 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("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") arrows(rep(0,p),rep(0,p),0.92*Atot[,axespar[1]],0.92*Atot[,axespar[2]],col=colpar+5) abline(v=0) abline(h=0) text(Atot[, axespar], dimnames(Atot)[[1]],col=colpar+5) if(vsup==1){ arrows(rep(0,nrow(As)),rep(0,nrow(As)),0.92*As[,axespar[1]],0.92*As[,axespar[2]],col=colpar+4) text(As[, axespar,drop=F], dimnames(Xc)[[2]],col=colpar+4) } } } else { cat("cercle des correlations (o/n) ?\n") plta <- scan("", character(), 1) if(length(plta) == 0) plta<-"n" if((plta == "N") | (plta == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltah <- scan("", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltav <- scan("", 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("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), "%)")) 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+5) text(Atot[, axespar], dimnames(Atot)[[1]],col=colpar+5) if(vsup==1){ arrows(rep(0,nrow(As)),rep(0,nrow(As)),0.92*As[,axespar[1]],0.92*As[,axespar[2]],col=colpar+4) text(As[, axespar,drop=F], dimnames(Xc)[[2]],col=colpar+4) } } } } # fin des graphiques } invisible( return(nomfichX, Xini, X, Q, D, inertiaX, k, valp, Fa, A, C, CTRus, COSus, CTRva, COSva, As, 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="lavender",cexpar=0.7) { # # 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 "cornsilk" ou "lightgray" utiliser aussi colours()[624] # # 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)) 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) cat(paste("D2 d'independance =", format(chi2 ),", d.d.l. =", ddl, ", Chi2 critique =", format( qchisq(0.95, ddl)),"(0.05)"), "\n") cat("_______________________________________________________________\n" ) cat("diagramme des valeurs propres (o ou n) ?\n") plth <- scan("", character(), 1) if(length(plth) == 0) plth <- "n" if(plth == "o" || plth == "O") { 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){ par(bg=bgpar) palette(rainbow(12)) par(mfrow = c(1, 1)) barplot(valp, space = 2, names = format(1:J), col = 1:max(I,J) ) title("valeurs propres") palette(rainbow(6)) } } } # 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("", 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("", 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("", character(), 1) if(length(pltc) == 0)pltc<-"n" if((pltc == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltch <- scan("", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltcv <- scan("", 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("", 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("", 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("",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("", character(), 1) if(length(plta) == 0)plta<-"n" if((plta == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltah <- scan("", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltav <- scan("", 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("", 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("",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("", character(), 1) if(length(plta) == 0)plta<-"n" if((plta == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltah <- scan("", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltav <- scan("", 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("", 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("", 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("",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("", 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("",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("", 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("",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("", 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("", numeric(), 1) cat("axe horizontal (<=", k, ") ?\n") pltah <- scan("", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltav <- scan("", 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 colories selon les modalites d'une variable qualitative? (o/n)\n") plto <- scan("", 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("", 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("",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("", 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("",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(nomfichN, Nini, N, DI, DJ, inertia, k, valp, A, C, CTRus, COSus, CTRva, COSva,coordlignesup,coordcolsup) if((!missing(colsup))& missing(lignesup))return(nomfichN, Nini, N, DI, DJ, inertia, k, valp, A, C, CTRus, COSus, CTRva, COSva,coordcolsup) if((!missing(lignesup))& missing(colsup)) return(nomfichN, Nini, N, DI, DJ, inertia, k, valp, A, C, CTRus, COSus, CTRva, COSva,coordlignesup) if(missing(colsup)& missing(lignesup)) return(nomfichN, Nini, N, DI, DJ, inertia, k, valp, A, C, CTRus, COSus, CTRva, 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("", 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("", 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("", 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("", character(), 1) if((length(pltc) == 0) | (pltc == "n")) break else { cat("axe horizontal (<=", k, ") ?\n") pltch <- scan("", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltcv <- scan("", 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("", 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("", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltcv <- scan("", numeric(), 1) cat("numero de la var. qual. (<=", q, ") ?\n") pltvqual <- scan("", 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("", 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("", 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("", 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("", numeric(), 1) cat("axe vertical (<=", k, ") ?\n") pltcv <- scan("", numeric(), 1) cat("numero de la var. en Z (<=", k, ") ?\n") pltcz <- scan("", 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(X,nbmod,U,B,labels) } #################################################### affecte<-function(X,g,Y,imin=1,imax=nrow(Y),D=1) { # entrees: # X matrice dont les lignes sont les coordonnees de points # g vecteur des modalites des groupes # Y matrice des vecteurs lignes des coordonnees des points a affecter; # 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) # 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) 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(missing(Y)) { for(i in 1:nrow(X)) tableau[groupe[i],g[i]]<-tableau[groupe[i],g[i]]+1 dimnames(tableau)<-list(paste("gaf",format(1:ng)),paste("gr",format(1:ng))) return(A,V,groupe,tableau) } else return(A,V,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) { # découpage MultiVariables de variables quantitatives # 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 #SORTIES #X matrice d'entiers des modalites #U matrice du codage disjonctif complet #B tableau de Burt #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 # # X<-as.matrix(X) n<-nrow(X) pX<-ncol(X) if(!missing(qual)){ if(qual>1){ Xquant<-X[,1:(qual-1),drop=F] } 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]) } } } else {nbmod<-rep(0,p) for(i in 1:p) {breaks[[i]]<-c(breaks[[i]],max(Xquant[,i])) nbmod[i]<-length(breaks[[i]]) } } 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("", 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("",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 if(is.matrix(Xquant)) return(X,U,B,nbmod,labels,breaks) else return(U,B,nbmod,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()) } #################################################### 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("",numeric(),1) cat("Nombre de lignes ?\n") lign<-scan("",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("",numeric(),1) cat("Nombre de lignes ?\n") lign<-scan("",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(P,margeL,margeC,ProfilsL,ProfilsC) } ####################################################