|
Solutions pour la séance numéro 3 (énoncés)
|
Retour à la page principale de
(gH)
Programmer en R
Session de formation, décembre 2013
gilles.hunault "at" univ-angers.fr
Solutions pour la séance numéro 3 (énoncés)
La fonction read.table() contient de nombreux paramètres qui couvrent la plupart des cas classiques. En particulier pour les noms de colonne en ligne 1, il y a header=TRUE. Et pour les noms de ligne, on utilise row.names(). Enfin, avec encoding="latin1" on peut lire les caractères accentués. Avec un peu d'habitude, on finit par savoir s'en servir :
read.table() usage: ------------------- read.table(file, header = FALSE, sep = "", quote = "\"'", dec = ".", row.names, col.names, as.is = !stringsAsFactors, na.strings = "NA", colClasses = NA, nrows = -1, skip = 0, check.names = TRUE, fill = !blank.lines.skip, strip.white = FALSE, blank.lines.skip = TRUE, comment.char = "#", allowEscapes = FALSE, flush = FALSE, stringsAsFactors = default.stringsAsFactors(), fileEncoding = "", encoding = "unknown", text)Pour retirer la dernière ligne, une fois les données lues, on utilise l'index négatif de cette ligne. La fonction row.names() prend un vecteur de noms uniques pour en faire des noms de lignes. Les constructions classiques de noms de ligne utilisent donc paste(), substr() et sprintf().
# si la ligne 1 contient le nom des lignes, la lecture simple ne suffit pas cats("une lecture ratée") data <- read.table("derniere.txt") print(data) # donc il faut ajouter header=TRUE cats("une lecture réussie") data <- read.table("derniere.txt",header=TRUE,row.names=1,encoding="latin1") data <- data[ -nrow(data), ] print(data) # normalisation des noms de ligne cats("aménagement des noms de ligne (1)") row.names(data) <- paste("Lig_",sprintf("%03d",as.numeric(row.names(data))),sep="") print(data) # une autre définition des noms de ligne cats("aménagement des noms de ligne (2)") debutGenre <- substr(data[,"genre"],1,3) debutEspece <- substr(data[,"espèce"],1,3) row.names(data) <- paste(debutGenre,"_",debutEspece,"_", sprintf("%03d",1:nrow(data)),sep="") print(data)une lecture ratée ================= V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 1 num accession length reign pfam cdd foldindex pi mw gravy genre esp\xe8ce 2 35 ABS44867 200 Viridiplantae N/A N/A -0.0835869 6.1399999 20295.8008 -1.0180905 Oryza sativa 3 822 A2211386C 153 Viridiplantae N/A N/A -0.1659534 9.3100004 15909.3203 -1.2745098 Pseudotsuga menziesii 4 964 ABB55135 132 Viridiplantae N/A N/A -0.2853720 9.4399996 14103.7598 -1.5803030 Pinus yunnanensis 5 965 AAW59164 143 Viridiplantae N/A N/A -0.1855139 8.6999998 15209.6104 -1.3573427 Pinus taeda 6 966 AAW59168 144 Viridiplantae N/A N/A -0.2060405 9.3100004 15422.8604 -1.4013889 Pinus taeda 7 980 AAW59174 143 Viridiplantae N/A N/A -0.1803197 8.6999998 15237.6699 -1.3405594 Pinus taeda 8 981 AAW59176 143 Viridiplantae N/A N/A -0.2015961 9.3100004 15294.7100 -1.3867133 Pinus taeda 9 total 7 7 7 0 0 7 7 7 7 7 7 une lecture réussie =================== accession length reign pfam cdd foldindex pi mw gravy genre espèce 35 ABS44867 200 Viridiplantae N/A N/A -0.0835869 6.14 20295.80 -1.018091 Oryza sativa 822 A2211386C 153 Viridiplantae N/A N/A -0.1659534 9.31 15909.32 -1.274510 Pseudotsuga menziesii 964 ABB55135 132 Viridiplantae N/A N/A -0.2853720 9.44 14103.76 -1.580303 Pinus yunnanensis 965 AAW59164 143 Viridiplantae N/A N/A -0.1855139 8.70 15209.61 -1.357343 Pinus taeda 966 AAW59168 144 Viridiplantae N/A N/A -0.2060405 9.31 15422.86 -1.401389 Pinus taeda 980 AAW59174 143 Viridiplantae N/A N/A -0.1803197 8.70 15237.67 -1.340559 Pinus taeda 981 AAW59176 143 Viridiplantae N/A N/A -0.2015961 9.31 15294.71 -1.386713 Pinus taeda aménagement des noms de ligne (1) ================================= accession length reign pfam cdd foldindex pi mw gravy genre espèce Lig_035 ABS44867 200 Viridiplantae N/A N/A -0.0835869 6.14 20295.80 -1.018091 Oryza sativa Lig_822 A2211386C 153 Viridiplantae N/A N/A -0.1659534 9.31 15909.32 -1.274510 Pseudotsuga menziesii Lig_964 ABB55135 132 Viridiplantae N/A N/A -0.2853720 9.44 14103.76 -1.580303 Pinus yunnanensis Lig_965 AAW59164 143 Viridiplantae N/A N/A -0.1855139 8.70 15209.61 -1.357343 Pinus taeda Lig_966 AAW59168 144 Viridiplantae N/A N/A -0.2060405 9.31 15422.86 -1.401389 Pinus taeda Lig_980 AAW59174 143 Viridiplantae N/A N/A -0.1803197 8.70 15237.67 -1.340559 Pinus taeda Lig_981 AAW59176 143 Viridiplantae N/A N/A -0.2015961 9.31 15294.71 -1.386713 Pinus taeda aménagement des noms de ligne (2) ================================= accession length reign pfam cdd foldindex pi mw gravy genre espèce Ory_sat_001 ABS44867 200 Viridiplantae N/A N/A -0.0835869 6.14 20295.80 -1.018091 Oryza sativa Pse_men_002 A2211386C 153 Viridiplantae N/A N/A -0.1659534 9.31 15909.32 -1.274510 Pseudotsuga menziesii Pin_yun_003 ABB55135 132 Viridiplantae N/A N/A -0.2853720 9.44 14103.76 -1.580303 Pinus yunnanensis Pin_tae_004 AAW59164 143 Viridiplantae N/A N/A -0.1855139 8.70 15209.61 -1.357343 Pinus taeda Pin_tae_005 AAW59168 144 Viridiplantae N/A N/A -0.2060405 9.31 15422.86 -1.401389 Pinus taeda Pin_tae_006 AAW59174 143 Viridiplantae N/A N/A -0.1803197 8.70 15237.67 -1.340559 Pinus taeda Pin_tae_007 AAW59176 143 Viridiplantae N/A N/A -0.2015961 9.31 15294.71 -1.386713 Pinus taedaPour avoir une numérotation interne au couple genre/espèce, le plus simple est de trier par couple et ensuite de numéroter dans une boucle :
# tri genre espèce pour exercice 1 séance 3 data <- read.table("derniere.txt",header=TRUE,row.names=1,encoding="latin1",stringsAsFactors=FALSE) data <- data[ -nrow(data), ] dao <- data[order(da$genre,da$espèce),] dao$genre <- reformate(dao$genre) dao$espèce <- reformate(dao$espèce) cpge <- function(x) { nbe <- length(x) # voici ce qu'il ne faut pas faire : npc <- rep(1,nbe) # écrire une fonction sans commentaire for (ide in (2:nbe)) { # avec des noms de variable non explicités if (x[ide] ==x[ide-1]) { # npc[ide] <- npc[ide-1]+1 # bonjour la lisibilité ! } # fin si } # fin pour return(npc) } # fin de fonction debutGenre <- substr(dao[,"genre"],1,3) debutEspece <- substr(dao[,"espèce"],1,3) numge <- cpge(paste(dao$genre,dao$espèce)) nomslig <- paste(debutGenre,"_",debutEspece,"_", sprintf("%03d",numge),sep="") row.names(dao) <- paste(debutGenre,"_",debutEspece,"_", sprintf("%03d",numge),sep="") print(dao)accession length reign pfam cdd foldindex pi mw gravy genre espèce Ory_sat_001 ABS44867 200 Viridiplantae N/A N/A -0.0835869 6.14 20295.80 -1.018091 Oryza _____ sativa ____ Pin_tae_001 AAW59164 143 Viridiplantae N/A N/A -0.1855139 8.70 15209.61 -1.357343 Pinus _____ taeda _____ Pin_tae_002 AAW59168 144 Viridiplantae N/A N/A -0.2060405 9.31 15422.86 -1.401389 Pinus _____ taeda _____ Pin_tae_003 AAW59174 143 Viridiplantae N/A N/A -0.1803197 8.70 15237.67 -1.340559 Pinus _____ taeda _____ Pin_tae_004 AAW59176 143 Viridiplantae N/A N/A -0.2015961 9.31 15294.71 -1.386713 Pinus _____ taeda _____ Pin_yun_001 ABB55135 132 Viridiplantae N/A N/A -0.2853720 9.44 14103.76 -1.580303 Pinus _____ yunnanensis Pse_men_001 A2211386C 153 Viridiplantae N/A N/A -0.1659534 9.31 15909.32 -1.274510 Pseudotsuga menziesii _Avec les paramètres skip, comment.char, blank.lines.skip et nrows de la fonction read.table(), lire partiellement un fichier avec des lignes vides et des commentaires est un jeu d'enfant :
data <- read.table( file="videcmt.txt", skip=5, comment.char="#", blank.lines.skip=TRUE, row.names=1, nrows=8 ) # fin de read.table colnames(data) <- c("Duree","Bolus","Serie") print(data)Duree Bolus Serie A35 617 552 1 A37 417 652 1 A38 410 565 1A A39 390 700 1A A40 500 540 1A A60 500 540 1B A61 480 408 1B A92 480 408 2Avec les indices dans les crochets et la séquence avec le symbole deux-points, on peut afficher des parties de dataframe. Si on utilise ncol() et nrow(), on peut trouver les indices des dernières lignes et colonnes, mais les fonctions head() et tail() sont à l'usage plus pratiques. Pour transposer, il faut utiliser la fonction t().
# lecture des données dans le "mauvais sens" cats("lecture") data <- read.table("transp.txt",header=TRUE,row.names=1) print(data) # transposition cats(" on transpose") data <- t(data) print(data) # premières/dernières lignes/colonnes cats("affichages") nbLig <- nrow(data) nbCol <- ncol(data) colnames(data) <- paste(colnames(data),"_V",1:nbCol,sep="") print( data[ (1:3), ] ) # ou head(data,n=3) cat("...\n") print( data[ ((nbLig-3):nbLig), (1:3)] ) # ou tail(data,n=4) cat("\n") print( data[ (3:5) , (1:3) ] ) print( data[ (5:7) , ((nbCol-2):nbCol) ] )lecture ======= M001 M002 M003 M004 M005 M006 M007 M008 M009 M010 SEXE 1 0 1 1 0 1 1 1 0 1 AGE 62 60 31 27 22 70 19 53 62 63 PROF 1 9 9 8 8 4 8 6 16 16 ETUD 2 3 4 4 4 1 4 2 4 0 REGI 2 4 4 1 1 1 4 2 2 1 USAG 3 1 1 1 2 1 2 3 2 0 on transpose ============= SEXE AGE PROF ETUD REGI USAG M001 1 62 1 2 2 3 M002 0 60 9 3 4 1 M003 1 31 9 4 4 1 M004 1 27 8 4 1 1 M005 0 22 8 4 1 2 M006 1 70 4 1 1 1 M007 1 19 8 4 4 2 M008 1 53 6 2 2 3 M009 0 62 16 4 2 2 M010 1 63 16 0 1 0 affichages ========== SEXE_V1 AGE_V2 PROF_V3 ETUD_V4 REGI_V5 USAG_V6 M001 1 62 1 2 2 3 M002 0 60 9 3 4 1 M003 1 31 9 4 4 1 ... SEXE_V1 AGE_V2 PROF_V3 M007 1 19 8 M008 1 53 6 M009 0 62 16 M010 1 63 16 SEXE_V1 AGE_V2 PROF_V3 M003 1 31 9 M004 1 27 8 M005 0 22 8 ETUD_V4 REGI_V5 USAG_V6 M005 4 1 2 M006 1 1 1 M007 4 4 2Pour filtrer des noms de ligne, la fonction regexpr() est tout-à-fait adaptée :
cats("toutes les données") data <- read.table("pas_2_3.txt",head=TRUE,row.names=1) print(data) cats("filtrage") flt <- regexpr("_[23]$",row.names(data))>0 print(data[ !flt ,])toutes les données ================== Duree Bolus Serie A35 617 552 1 A37 417 652 1 A38_1 410 565 1A A39 390 700 1A A40_2 500 540 1A A60_2 500 540 1B A61 480 408 1B A92_3 480 408 2 filtrage ======== Duree Bolus Serie A35 617 552 1 A37 417 652 1 A38_1 410 565 1A A39 390 700 1A A61 480 408 1BPour des affichages sélectifs de lignes, les conditions logiques avec & et | suffisent souvent :
elf <- lit.dar.wd("elf.dar") cats("femmes de moins de 30 ans (début)") flt <- (elf$SEXE==1) & (elf$AGE<30) head( elf[ flt , c("SEXE","AGE") ] ) cats("hommes ou personnes de plus de 60 ans (début)") flt <- (elf$SEXE==0) | (elf$AGE>60) head( elf[ flt , c("SEXE","AGE") ] )femmes de moins de 30 ans (début) ================================= SEXE AGE M004 1 27 M007 1 19 M012 1 11 M019 1 21 M020 1 23 M025 1 14 hommes ou personnes de plus de 60 ans (début) ============================================= SEXE AGE M001 1 62 M002 0 60 M005 0 22 M006 1 70 M009 0 62 M010 1 63Pour accéder à une colonne d'un dataframe, on peut utiliser son nom avec $, son numéro de colonne ou son nom comme indice. Il est conseillé d'utiliser le nom de colonne plutôt que le numéro parce que c'est plus sûr. Pour supprimer, on peut utiliser l'indice négatif mais la notation en liste avec dollar et NULL est plus concise.
elf <- lit.dar.wd("elf.dar") mini <- head(elf) cats("trois moyens d'accéder à la colonne AGE") v1 <- mini[ ,2 ] v2 <- mini[ ,"AGE" ] v3 <- mini$AGE print( cbind( v1, v2, v3 ) ) cats("ce sont bien les mêmes structures :") cat(" v1=v2 ? ",identical(v1,v2),"\n") cat(" v1=v3 ? ",identical(v1,v3),"\n") cat(" v2=v2 ? ",identical(v2,v3),"\n") cats("supprimons les lignes 3 à 5") mini2 <- mini[ -(3:5), ] print(mini2) # pour supprimer la colonne AGE : # solution 1 mini <- mini[, -2] # solution 2 mini <- mini[, -which("AGE"==names(mini)) ] # solution 3 mini$AGE <- NULLtrois moyens d'accéder à la colonne AGE ======================================= v1 v2 v3 [1,] 62 62 62 [2,] 60 60 60 [3,] 31 31 31 [4,] 27 27 27 [5,] 22 22 22 [6,] 70 70 70 ce sont bien les mêmes structures : =================================== v1=v2 ? TRUE v1=v3 ? TRUE v2=v2 ? TRUE supprimons les lignes 3 à 5 =========================== SEXE AGE PROF ETUD REGI USAG M001 1 62 1 2 2 3 M002 0 60 9 3 4 1 M006 1 70 4 1 1 1Pour renommer occasionnellement une colonne, on peut utiliser son numéro et modifier le colnames() de la structure. Pour du renommage "massif", il vaut mieux passer par la fonction rename.pattern() du package epicalc.
> # les données > > elf <- lit.dar.wd("elf.dar") > mini <- head(elf) > print( mini ) SEXE AGE PROF ETUD REGI USAG M001 1 62 1 2 2 3 M002 0 60 9 3 4 1 M003 1 31 9 4 4 1 M004 1 27 8 4 1 1 M005 0 22 8 4 1 2 M006 1 70 4 1 1 1 > # on veut renommer la colonne "PROF" en "csp" > # solution (1) : on connait le numéro > > print(cbind(colnames(mini))) [,1] [1,] "SEXE" [2,] "AGE" [3,] "PROF" [4,] "ETUD" [5,] "REGI" [6,] "USAG" > colnames(mini)[3] <- "csp" > print( mini ) SEXE AGE csp ETUD REGI USAG M001 1 62 1 2 2 3 M002 0 60 9 3 4 1 M003 1 31 9 4 4 1 M004 1 27 8 4 1 1 M005 0 22 8 4 1 2 M006 1 70 4 1 1 1 > # on veut renommer la colonne "csp" en "prof" > # solution (2) : on ne connait pas le numéro > > avant <- "csp" > apres <- "prof" > noms <- colnames(mini) > colnames(mini)[ which(noms==avant) ] <- apres > print( mini ) SEXE AGE prof ETUD REGI USAG M001 1 62 1 2 2 3 M002 0 60 9 3 4 1 M003 1 31 9 4 4 1 M004 1 27 8 4 1 1 M005 0 22 8 4 1 2 M006 1 70 4 1 1 1 > # on veut renommer la colonne "prof" en "csp" > # solution (3) : rename de epicalc > > library(epicalc) Loading required package: foreign Loading required package: survival Loading required package: splines Loading required package: MASS Loading required package: nnet > rename.pattern(x1="prof",x2="csp",dataFrame=mini) Note the following change(s) in variable name(s): Old var names New var names "prof" "csp" > print( mini ) SEXE AGE csp ETUD REGI USAG M001 1 62 1 2 2 3 M002 0 60 9 3 4 1 M003 1 31 9 4 4 1 M004 1 27 8 4 1 1 M005 0 22 8 4 1 2 M006 1 70 4 1 1 1Le package XML et ses fonctions readHTMLTable() et readHTMLList() est prévu pour cela. Il faut sans doute juste utiliser des fonctions sur chaines comme iconv() et str_replace() afin d'obtenir des informations "propres".
## récupération d'un tableau XHTML dans R # chargement du package XML library(XML) # chargement du package stringr library(stringr) # pour la fonction str_replace() # url de la page Web à utiliser url <- "http://forge.info.univ-angers.fr/~gh/wstat/Programmation_R/progr.php?n=1&m=s" tabp <- readHTMLTable(url,which=3,stringsAsFactors=FALSE) print(tabp,right=FALSE) # un peu de "ménage" colnames(tabp) <- str_replace(colnames(tabp),"\\n","") # tout sauf les deux derniers caractères tsdc <- function(x) { substr(x,1,nchar(x)-2) } tabp[,2] <- as.numeric( tsdc(tabp[,2]) ) print(tabp,right=FALSE) ## récupération d'une liste ip <- "http://fr.wikipedia.org/wiki/Institut_Pasteur" cdr <- readHTMLList(ip,which=7) # conversion des accents iconv(cdr,from="UTF-8",to="LATIN1") print(cbind(cdr))> ## récupération d'un tableau XHTML dans R > > # chargement du package XML > > library(XML) > # chargement du package stringr > > library(stringr) # pour la fonction str_replace() > # url de la page Web à utiliser > > url <- "http://forge.info.univ-angers.fr/~gh/wstat/Programmation_R/progr.php?n=1&m=s" > tabp <- readHTMLTable(url,which=3,stringsAsFactors=FALSE) > print(tabp,right=FALSE) \nPackage \nNb_objets \nLien 1 base 1167Â lls_base.sor 2 gdata 64Â lls_gdata.sor 3 tools 97Â lls_tools.sor 4 utils 198Â lls_utils.sor > # un peu de "ménage" > > colnames(tabp) <- str_replace(colnames(tabp),"\\n","") > # tout sauf les deux derniers caractères > > tsdc <- function(x) { substr(x,1,nchar(x)-2) } > tabp[,2] <- as.numeric( tsdc(tabp[,2]) ) > print(tabp,right=FALSE) Package Nb_objets Lien 1 base 1167 lls_base.sor 2 gdata 64 lls_gdata.sor 3 tools 97 lls_tools.sor 4 utils 198 lls_utils.sor > ## récupération d'une liste > > ip <- "http://fr.wikipedia.org/wiki/Institut_Pasteur" > cdr <- readHTMLList(ip,which=7) > # conversion des accents > > iconv(cdr,from="UTF-8",to="LATIN1") [1] "Biologie cellulaire et infection" "Biologie du développement et cellules souches" [3] "Biologie structurale et chimie" "Génomes et génétique" [5] "Immunologie" "Infection et épidémiologie" [7] "Microbiologie" "Neuroscience" [9] "Parasitologie et mycologie" "Virologie" > print(cbind(cdr)) cdr [1,] "Biologie cellulaire et infection" [2,] "Biologie du développement et cellules souches" [3,] "Biologie structurale et chimie" [4,] "Génomes et génétique" [5,] "Immunologie" [6,] "Infection et épidémiologie" [7,] "Microbiologie" [8,] "Neuroscience" [9,] "Parasitologie et mycologie" [10,] "Virologie"Si ce genre de récupération vous intéresse, il y a d'autres exemples dans notre introduction non élémentaire à R, exercice 5 du cours 2.
Le code suivant est presque suffisant pour lire un entier ou des entiers au clavier :
litEntier <- function() { repeat { lu <- readline("litEntier> donner un nombre entier : ") lulu <<- lu # le doublement de < n'est pas une erreur try(n <- as.integer(lu)) if (is.na(n)) { cat("litEntier> hélas, ce n'est pas un entier !\n") } else { break } # fin de si } # fin de répéter cat("Merci. Votre nombre est ",n," et son carré est ",n*n,"\n") return( n) } # fin de fonction litEntier # ------------------------------------------------------------------------ litEntierAvecRegexpr <- function() { nbLect <- 0 repeat { if (nbLect==0) { demande <- "litEntierAvecRegexpr> donner un nombre entier : " } else { demande <- "litEntierAvecRegexpr> redonner un nombre entier : " } # fin si lu <- readline(demande) nbLect <- nbLect + 1 if ((regexpr("^[0-9]+$",lu,perl=TRUE)==1)) { # attention : "[0-9]+" n'est pas suffisant break } # fin de si cat("litEntierAvecRegexpr> hélas, ce n'est pas un entier !\n") } # fin de répéter n <- as.numeric(lu) cat("Merci. Votre nombre est ",n," et son carré est ",n*n,"\n") return( n) } # fin de fonction litEntierAvecRegexpr # ------------------------------------------------------------------------ litEntiers <- function() { repeat { lus <- readline("donner des nombre entiers séparés par des espaces : ") if ((regexpr("^[0-9 ]+$",lus,perl=TRUE)==1)) { break } # fin de si cat("hélas, ce ne sont pas des entiers séparés par des espaces !\n") } # fin de répéter lulu <<- lus # le doublement de < n'est pas une erreur ldm <- unlist( strsplit(x=lus,split="\\s+")) vect <- as.vector(sapply(X=ldm,FUN=as.integer)) cat("Merci. Vous avez entré ",length(vect)," entier(s).\n") return(vect) } # fin de fonction litEntiersdémo de litEntier ================= litEntier> donner un nombre entier : oui litEntier> hélas, ce n'est pas un entier ! litEntier> donner un nombre entier : 17 25 litEntier> hélas, ce n'est pas un entier ! litEntier> donner un nombre entier : -6.1 Merci. Votre nombre est -6 et son carré est 36 démo de litEntierAvecRegexpr ============================ litEntierAvecRegexpr> donner un nombre entier : oui litEntierAvecRegexpr> hélas, ce n'est pas un entier ! litEntierAvecRegexpr> redonner un nombre entier : 17 25 litEntierAvecRegexpr> hélas, ce n'est pas un entier ! litEntierAvecRegexpr> redonner un nombre entier : 6.1 litEntierAvecRegexpr> hélas, ce n'est pas un entier ! litEntierAvecRegexpr> redonner un nombre entier : -3 litEntierAvecRegexpr> hélas, ce n'est pas un entier ! litEntierAvecRegexpr> redonner un nombre entier : 3 Merci. Votre nombre est 3 et son carré est 9 démo de litEntiers ================== donner des nombre entiers séparés par des espaces : 5 12 17 Merci. Vous avez entré 3 entier(s). [1] 5 12 17 Messages d'avis : 1: In doTryCatch(return(expr), name, parentenv, handler) : NAs introduits lors de la conversion automatique 2: In doTryCatch(return(expr), name, parentenv, handler) : NAs introduits lors de la conversion automatique 3: In lapply(X = X, FUN = FUN, ...) : NAs introduits lors de la conversion automatiqueMais ce code n'est pas correct, ni efficace : passer par as.integer est maladroit, car pour 6.1 cela ne renvoie pas d'erreur. Et malheureusement, comme le montre l'affichage précédent, la première fonction nommée litEntier() génère des messages d'avertissement (et non pas des messages d'erreurs). Donc l'aide sur try() qui propose le paramètre silent=TRUE n'est pas adapté. Comme les "warnings" sont gérés par les options (chercher warn), le code suivant est plus "propre" :
litEntier <- function() { options(warn=-1) # mieux : sovWarn <- as.integer( options("warn") ) # puis options(warn=-1) repeat { lu <- readline("litEntier> donner un nombre entier : ") lulu <<- lu try(n <- as.integer(lu)) if (is.na(n)) { cat("litEntier> hélas, ce n'est pas un entier !\n") } else { break } # fin de si } # fin de répéter options(warn=0) # mieux : options(warn=sovWarn) cat("Merci. Votre nombre est ",n," et son carré est ",n*n,"\n") return( n) } # fin de fonction litEntierDe même, la fonction nommée litEntierAvecRegexpr() ne sait pas reconnaitre un entier signé. Rajouter juste le signe moins et donc utiliser ^[-0-9]$ ne suffit pas : +3 ne serait pas reconnu. Nous vous laissons comprendre la solution suivante.
litEntierAvecRegexprEtSigne <- function() { nbLect <- 0 repeat { if (nbLect==0) { demande <- "litEntierAvecRegexprEtSigne> donner un nombre entier : " } else { demande <- "litEntierAvecRegexprEtSigne> redonner un nombre entier : " } # fin si lu <- readline(demande) nbLect <- nbLect + 1 if ((regexpr("^[+-]?[0-9]+$",lu,perl=TRUE)==1)) { # attention : "[0-9]+" n'est pas suffisant break } # fin de si cat("litEntierAvecRegexprEtSigne> hélas, ce n'est pas un entier !\n") } # fin de répéter n <- as.numeric(lu) cat("Merci. Votre nombre est ",n," et son carré est ",n*n,"\n") return( n) } # fin de fonction litEntierAvecRegexprEtSigne # ------------------------------------------------------------------------Il s'agit d'une généralisation de l'exercice précédent. Pour des nombres réels éventuellement signés, l'expression régulière est un peu plus compliquée. Pour des fractions, une fois passé le cap de la validité de l'expression, il faut convertir chaque fraction en un nombre, par exemple avec notre fonction convFrac().
litNombres <- function(chaine) { library(stringr) if (missing(chaine)) { cat("paramètre manquant, la syntaxe est litNombres(chaine) \n"); cat("par exemple : litNombres(\"1 3 12.8 -7 3 \")\n\n") return() } # fin de si if (chaine=="") { return( numeric(0) ) } # fin de si if (regexpr("^( *[+-]?[0-9]*(\\.|,)?[0-9]* *)*$",chaine,perl=TRUE)==1) { cat("ok\n") chaine <- str_replace_all(chaine,",",".") vectn <- unlist(strsplit(x=ghtrim(chaine),split="\\s+",perl=TRUE)) return(as.numeric(vectn)) } else { cat("non\n") return( c() ) } # fin si } # fin de fonction litNombres # ---------------------------------------- convFrac <- function(chaine) { chaine <- ghtrim(chaine) # supprime les espaces de début et fin if (1!=regexpr("[0-9]+/[0-9]+",chaine)) { return(NA) } else { parts <- as.numeric(unlist(strsplit(chaine,"/"))) return(parts[1]/parts[2]) } # fin si } # fin de fonction convFrac # ---------------------------------------- litFractions <- function(chaine) { vectFract <- unlist(strsplit(x=ghtrim(chaine),split=" +",perl=TRUE)) cat("on est passé de : ") str(chaine) cat("à : ") str(vectFract) cat("\n") return(as.vector(sapply(vectFract,convFrac))) } # fin de fonction litFractions # ---------------------------------------- cats("exemples d'utilisations :" ) v <- litNombres("1 3 12.8 -7 3 ") print( str(v) ) fr <- "3/4" cat("la fraction ",fr," vaut ",convFrac(fr),"\n") (frs <- litFractions(" 4/3 2/3 1/5") )exemples d'utilisations : ========================= > v <- litNombres("1 3 12.8 -7 3 ") ok > print( str(v) ) num [1:5] 1 3 12.8 -7 3 NULL > w <- litNombres("1 3 12.8 -7,5 3 ") ok > print(w) [1] 1.0 3.0 12.8 -7.5 3.0 > fr <- "3/4" > cat("la fraction ",fr," vaut ",convFrac(fr),"\n") la fraction 3/4 vaut 0.75 > (frs <- litFractions(" 4/3 2/3 1/5") ) on est passé de : chr " 4/3 2/3 1/5" à : chr [1:3] "4/3" "2/3" "1/5" [1] 1.3333333 0.6666667 0.2000000Contrairement à <-, l'opérateur <<- réalise une affectation dans l'environnement global. Donc si on écrit n <- 2 dans une fonction, n est une variable locale et n'existe pas en dehors (ou, si n existait avant, elle est masquée dans la fonction par la variable n locale, mais non modifiée). Par contre, si on écrit n <<- 2 dans une fonction, n devient une variable globale et remplace éventuellement celle qui existait déjà :
queVautn <- function() { if (exists("n")) { cat("n vaut ",n,"\n") } else { cat("n n'existe pas.\n") } # fin si } # fin de fonction queVautn demoLocale <- function() { n <- 3 cat(" dans la fonction, n vaut ",n,"\n") } # fin de fonction demoLocale # démonstration n <- 1 queVautn() demoLocale() queVautn() rm(n) demoLocale() queVautn()n vaut 1 dans la fonction, n vaut 3 n vaut 1 dans la fonction, n vaut 3 n n'existe pas.Voici une façon de s'en servir : vous êtes en train de mettre au point une fonction et vous voulez tester en ligne de commande des instructions sur des variables créées par la fonction. Utiliser <<- vous permet d'accéder à ces variables en dehors de la fonction...
Notre fonction analexies() de notre fichier statgh.r fait le travail demandé en ce qui concerne les dictionnaires :
# le fichier statgh.r contient la fonction analexies source("statgh.r") ## autres possibilités (changer éventuellement les chemins d'accès) : # source("/home/moi/stats/statgh.r") # source("http://forge.info.univ-angers.fr/~gh/statgh.r") # source("K:/Stat_ad/statgh.r") adt <- analexies( lit.texte("animaux.txt") ) print(names(adt)) print(adt)19 mots en tout dont 14 mots distincts. chaque mot est donc répété 1.36 fois en moyenne. [1] "nmots" "hmots" "nmotsd" "tmots" "toccs" $nmots [1] 19 $hmots <hash> containing 14 key-value pair(s). aboie : 1 alors : 1 animaux : 1 avec : 1 chat : 2 chien : 2 des : 1 et : 1 le : 4 miaule : 1 pattes : 1 quatre : 1 que : 1 sont : 1 $nmotsd [1] 14 $tmots occ aboie 1 alors 1 animaux 1 avec 1 chat 2 chien 2 des 1 et 1 le 4 miaule 1 pattes 1 quatre 1 que 1 sont 1 $toccs occ le 4 chien 2 chat 2 sont 1 que 1 quatre 1 pattes 1 miaule 1 et 1 des 1 avec 1 animaux 1 alors 1 aboie 1Pour une matrice d'occurences de mots dans des textes, on peut utiliser notre fonction cpt() dont voici un exemple d'utilisation
cpt("bete.txt candide.txt","maison temps train")# au bout d'environ 20 secondes : fichier 1 / 2 : bete.txt 23029 mots en tout dont 4238 mots distincts. chaque mot est donc répété 5.43 fois en moyenne. fichier 2 / 2 : candide.txt 1958 mots en tout dont 731 mots distincts. chaque mot est donc répété 2.68 fois en moyenne. maison temps train bete.txt 21 9 56 candide.txt 3 1 0Voici le code-source des fonctions et des exemples d'utilisation pour bien comprendre comment elles fonctionnnent.
prems renvoie tous les "premiers" éléments d'un vecteur, c'est-à-dire en fait tous les éléments sauf le dernier.
derns renvoie tous les "derniers" éléments d'un vecteur, c'est-à-dire tous les éléments sauf le premier.
affRec affiche récursivement la chaine de caractère qui représente formellement la somme des éléments d'un vecteur. Cette fonction affiche la somme formelle des premiers éléments, obtenue récursivement via prems(), puis le signe plus et enfin le dernier élément.
Contrairement à affRec, les fonctions somRecCrois et somRecFin calculent effectivement de façon récursive la somme des éléments d'un vecteur. somRecCrois effectue le calcul à partir du premier élément alors que somRecFin commence le calcul à partir du dernier élément.
prems <- function(vect) { # renvoie tous les éléments sauf le dernier return(vect[-length(vect)]) } # fin de fonction prems derns <- function(vect) { # renvoie tous les éléments sauf le premier return(vect[-1]) } # fin de fonction derns # exemple d'affichage récursif utilisant prem() affRec <- function(vect) { if(length(vect)==1) { return(vect) } else { return(paste(affRec(prems(vect))," + ",vect[length(vect)],sep="")) } # fin de si } # fin de fonction affRec somRecCrois <- function(vect,dbg=FALSE) { # somme récursive des éléments de vect à partir du premier if (dbg) { cat("appel de somRecCrois(",vect,")\n") } if(length(vect)==1) { return(vect) } else { return(vect[1]+somRecCrois(derns(vect),dbg)) } # fin de si } # fin de fonction somRecCrois somRecFin <- function(vect,dbg=FALSE) { # somme récursive des éléments de vect à partir du dernier if (dbg) { cat("appel de somRecFin(",vect,")\n") } if(length(vect)==1) { return(vect) } else { return(somRecFin(prems(vect),dbg)+vect[length(vect)]) } # fin de si } # fin de fonction somRecFin # vecteur de données et exemples v <- 10*(1:5) cats("Vecteur v") cat(v,"\n") cats("exemple d'affichage récursif ") cat(affRec(v),"\n") cats("Vecteur derns(v)") cat(derns(v),"\n") cats("Vecteur prems(v)") cat(prems(v),"\n") cats("Somme des éléments de v :") cat(sum(v),"\n") cats("Somme récursive par le début des éléments de v :") cat(somRecCrois(v),"\n") cats("explications") cat(somRecCrois(v,dbg=TRUE),"\n") cats("Somme récursive par la fin des éléments de v :") cat(somRecFin(v),"\n") cats("explications") cat(somRecFin(v,dbg=TRUE),"\n") # on ajoute des parenthèses pour y voir plus clair somRecCroisV2 <- function(vect) { if(length(vect)==1) { return(vect) } else { return(paste(vect[1]," + (",somRecCroisV2(derns(vect)),")",sep="")) } # fin de si } # fin de fonction somRecCroisV2 cats("détails 1") cat(somRecCroisV2(v),"\n") somRecFinV2 <- function(vect) { if(length(vect)==1) { return(vect) } else { return(paste("(",somRecFinV2(prems(vect)),") + ",vect[length(vect)],sep="")) } # fin de si } # fin de fonction somRecFinV2 cats("détails 2") cat(somRecFinV2(v),"\n") # tri récursif d'un vecteur : on trouve le minimum, on l'enlève # et on trie ce qui reste v <- round(runif(5,1,10),1) cats(" vecteur de départ : ") cat(v,"\n") triRec <- function(vect) { if (length(vect)==1) { return(vect) } else { mini <- min(vect) indMin <- which(vect==mini) return(c(mini,triRec(vect[-indMin]))) # possible, mais moins lisible : # } else { # mini <- min(vect) # return(c(mini,triRec(vect[-which(vect==mini)]))) } # fin de si } # fin de fonction triRec triRecV2 <- function(vect) { if (length(vect)==1) { return(vect) } else { mini <- min(vect) indMin <- which.min(vect) return(c(mini,triRec(vect[-indMin]))) } # fin de si } # fin de fonction triRecV2 w <- triRec(v) cats(" vecteur trié : ") cat(w,"\n") # démonstrations des tris récursifs v <- c(2,8,9,2,7,5,8) cat("\nAttention : le triRec de ",v,"est",triRec(v),"\n") cat("\nle triRecV2 de ",v,"est",triRecV2(v),"\n")Vecteur v ========= 10 20 30 40 50 exemple d'affichage récursif ============================= 10 + 20 + 30 + 40 + 50 Vecteur derns(v) ================ 20 30 40 50 Vecteur prems(v) ================ 10 20 30 40 Somme des éléments de v : ========================= 150 Somme récursive par le début des éléments de v : ================================================ 150 explications ============ appel de somRecCrois( 10 20 30 40 50 ) appel de somRecCrois( 20 30 40 50 ) appel de somRecCrois( 30 40 50 ) appel de somRecCrois( 40 50 ) appel de somRecCrois( 50 ) 150 Somme récursive par la fin des éléments de v : ============================================== 150 explications ============ appel de somRecFin( 10 20 30 40 50 ) appel de somRecFin( 10 20 30 40 ) appel de somRecFin( 10 20 30 ) appel de somRecFin( 10 20 ) appel de somRecFin( 10 ) 150 détails 1 ========= 10 + (20 + (30 + (40 + (50)))) détails 2 ========= ((((10) + 20) + 30) + 40) + 50 vecteur de départ : ===================== 9.5 5.1 6.9 4.5 8.3 vecteur trié : ================ 4.5 5.1 6.9 8.3 9.5 Attention : le triRec de 2 8 9 2 7 5 8 est 2 5 7 8 9 le triRecV2 de 2 8 9 2 7 5 8 est 2 2 5 7 8 9Comme on peut le voir, une fonction récursive a toujours un test d'arrêt (ou d'initialisation) et un appel à elle-même mais avec un paramètre de longueur inférieure. L'appel récursif peut être la première partie ou la dernière partie de la fonction.
Avec which(vect==minVec) on trouve (donc on peut supprimer) toutes les occurences du minimum alors que which.min(vect) ne fournit que la première occurence du minimum. Avec which(vect==minVec) on obtient donc un vecteur trié sans doublon, ce qui peut aussi se faire à l'aide des fonctions unique() et duplicated().
Il s'agit de notre fonction decritQT() de statgh.r, bien sûr. Elle a aussi un paramètre de langue, une valeur maximale en y pour l'histogramme des classes etc. Le tracé tige et feuilles et les graphiques sont gérés par l'appel de la fonction plotQTdet().
R est prévu pour s'exécuter en mode dit "BATCH". Donc un script de commandes pour le système d'exploitation suffit à lancer R, exécuter le programme et produire les fichiers de résultats. Toutefois, pour une meilleure portabilité, par exemple en cas de machines différentes (Windows, Linux, MacOs...) il est certainement préférable d'écrire un script dans un langage interactif comme perl ou rexx. Voici, à titre d'exemple notre fichier rstat.rex et des exemples d'utilisation
/* (gH) -- rstat.rex ; TimeStamp (unix) : 09 Janvier 2003 14:06 */ /* ceci exécute R avec un programme passé en paramètre. */ /* on peut préciser le nom du fichier de sortie */ parse arg nomprog nomsor if words(nomprog)=0 then do say say " syntaxe : rstat [ prog [fichier_sor] ] " say say " on lance R --no-save < $1 > $ 2 " say " rstat r lance R en interactif" say return end /* fin de si */ cmd = " R --no-save --encoding='latin1' " if (translate(nomprog)\="R") then cmd = cmd " < " nomprog if words(nomsor)>0 then do cmd = cmd " > " nomsor end /* fin de si */ say cmd cmd# appel de R en ligne de commande rstat r # utilisation d'un programme R nommé test.r en ligne de commandes # (afichage à l'écran) rstat test.r # utilisation d'un programme R nommé prog.r en ligne de commandes # (afichage dans le fichier prog.sor) rstat prog.r prog.sorSi on doit toujours exécuter disons la fonction allQT() sur un jeu de données, le plus simple est d'écrire un script R qui lit toujours le même fichier de données, puis d'adapter le script d'appel pour qu'il recopie le fichier des données de l'utilisateur avec le nom imposé avant d'exécuter ce script R...
Voici par exemple le programme qtslt.r à utiliser pour analyser le fichier qtslt.dar qui ne contient que des variables quantitatives :
########################################### # ## analyse du fichier qtslt.dar # ########################################### # 1. chargement des fonctions (gH) source("statgh.r",encoding="latin1") # 2. lecture des données qtslt <- lit.dar("qtslt.dar") # 3. appel de la fonction allQTdf allQTdf(qtslt)Pour utiliser ce programme et analyser vins.dar, voici ce qu'il faut faire sous Windows en ligne de commande :
copy vins.dar qtslt.dar "C:\Program Files\R-2.14\bin\i386\Rterm.exe" --no-save < D:\Cours\R\qtslt.r > vins_qtslt.txtEt sous Linux (ou MacOs) :
cp vins.dar qtslt.dar /usr/bin/R --no-save < ~/Cours/R/qtslt.r > vins_qtslt.txtPour utiliser ce programme et analyser iris.dar, voici ce qu'il faut faire sous Windows en ligne de commande :
copy iris.dar qtslt.dar "C:\Program Files\R-2.14\bin\i386\Rterm.exe" --no-save < D:\Cours\R\qtslt.r > iris_qtslt.txtEt sous Linux (ou MacOs) :
cp iris.dar qtslt.dar /usr/bin/R --no-save < ~/Cours/R/qtslt.r > iris_qtslt.txtCes commandes à écrire ne sont pas simples et il vaut mieux les éviter, sauf si on est un(e) professionnel(le) de l'informatique. Parmi les défauts qu'on peut leur reprocher, il y a le fait qu'il est difficile de se souvenir des chemins d'accès des logiciels. De plus, il n'est pas toujours évident de maitriser la redirection des entrées avec le caractère < et la redirection des sorties avec le caractère >. Enfin, les copies de noms de fichier "à la main" peuvent être parfois source d'erreurs.
Pour toutes ces raisons, le plus simple est d'écrire un programme qui s'en occupe à notre place. Voici une solution Rexx/Regina
/* (gH) -_- anaqt_rex.rex ; TimeStamp (unix) : 01 Mars 2012 vers 14:42 */ version = 1.0 parse arg nomf /* rappel de la syntaxe si aucun paramètre */ if words(nomf)=0 then do say say " anaqt_rex : analyse d'un fichier qui ne contient que des QT" say " syntaxe : anaqt_rex fichier[.dar] " say " exemples : anaqt_rex demo1.dar" say " anaqt_rex demo2 " say return end /* fin de si */ /* préparation des noms de fichier */ /* si le fichier se termine par ., on le retire */ if (substr(nomf,length(nomf),1)=".") then nomf = substr(nomf,1,length(nomf)-1) /* si le fichier ne se termine pas par .dar (à cause de la complétion via tab) */ /* on ajoute .dar au nom de fichier */ if (index(nomf,".dar")=0) then nomf = nomf||".dar" /* on vérifie que le fichier existe */ rf = fileut("checkFile^"nomf) if rf = -1 then do say " en conséquence, j'arrête le programme." return end /* fin de si */ /* on recopie le fichier de données en qtslt.dar*/ cmd = "cp " nomf " qtslt.dar" say cmd cmd /* on exécute le programme R avec une sortie dans le fichier standard */ cmd = "rstat qtslt.r qtslt.sor " say cmd cmd /* on recopie le fichier de résultats avec un nom plus sympathique */ bn = fileut("fileName^"nomf) fs = bn||".sor" cmd = "cp qtslt.sor " fs say cmd cmd /* affichage de fin */ say " vous pouvez utiliser " fsEt une solution Perl
# (gH) -_- anaqt_pl.pl ; TimeStamp (unix) : 01 Mars 2012 vers 15:41 # test du paramètre if ($ARGV[0] eq "") { print "\n" ; print " anaqt_pl : analyse d'un fichier qui ne contient que des QT \n" ; print " syntaxe : anaqt_pl fichier[.dar] \n" ; print " exemples : anaqt_pl demo1.dar \n" ; print " anaqt_pl demo2 \n" ; print "\n" ; exit(-1) ; } # fin de si $nomf = $ARGV[0] ; # préparation des noms de fichier # si le fichier se termine par ., on le retire if (substr($nomf,length($nomf)-1) eq ".") { $nomf = substr($nomf,0,length($nomf)-1) ; } ; # si le fichier ne se termine pas par .dar (à cause de la complétion via tab) # on ajoute .dar au nom de fichier if (index($nomf,".dar")==-1) { $nomf = $nomf.".dar" ; } ; # on vérifie que le fichier existe * open(FE,"<$nomf") or die("le fichier $nomf n'existe pas.\n\n") ; close(FE) ; # on recopie le fichier de données en qtslt.dar $cmd = "cp $nomf qtslt.dar" ; print "$cmd\n" ; system($cmd) ; # on exécute le programme R avec une sortie dans le fichier standard $cmd = "rstat qtslt.r qtslt.sor " ; print "$cmd\n" ; system($cmd) ; # on recopie le fichier de résultats avec un nom plus sympathique $bn = substr($nomf,0,rindex($nomf,".")) ; $fs = $bn.".sor" ; $cmd = "cp qtslt.sor $fs " ; print "$cmd\n" ; system($cmd) ; # affichage de fin print "\n vous pouvez utiliser $fs \n\n"Enfin, si on écrit un script exécutable qui exécute ces programmes, et qui est accessible par le PATH, l'utilisation est encore facilitée, sous Windows, comme sous Linux ou MacOs. Par exemple sous Linux, si on rend exécutable via chmod +x anaqt_rex le fichier suivant nommé anaqt_rex
regina ~gh/Bin/anaqt_rex.rex $*alors on peut utiliser directement la commande anaqt_rex. En voici la démonstration :
# démonstrations de l'utilisation du script Rexx os> anaqt_rex anaqt_rex : analyse d'un fichier qui ne contient que des QT syntaxe : anaqt_rex fichier[.dar] exemples : anaqt_rex demo1.dar anaqt_rex demo2 os> anaqt_rex pommes.dar -- ! -- 'Fichier' pommes.dar 'non vu dans le cd courant !' en conséquence, j'arrête le programme. os> anaqt_rex vins.dar cp vins.dar qtslt.dar rstat qtslt.r qtslt.sor cp qtslt.sor vins.sor vous pouvez utiliser vins.sor os> anaqt_rex iris cp iris.dar qtslt.dar rstat qtslt.r qtslt.sor cp qtslt.sor iris.sor vous pouvez utiliser iris.sorIl est clair qu'au vu de la démonstration du package hwriter que l'on pourra consulter ici, fournir les résultats en HTML est une bonne idée puisque les traitements de texte moderne (post 2010) savent ouvrir les fichiers HTML.
De façon générale, un navigateur permet de visualiser n'importe quel texte, permet de grossir ou de rétrécir la taille des caractères, c'est donc un bon candidat pour afficher des données importantes. C'est pourquoi au lieu de la fonction print() de R, nous préférons utiliser notre fonction hprint() qui affiche les structures dans le navigateur.
La solution consiste à compter le nombre de valeurs distinctes. Il y a en général peu de valeurs pour une QL, et souvent beaucoup plus pour une QT. Pour compter le nombre de valeurs distinctes du vecteur v, on peut utiliser length( unique(v) ). Voici donc la nouvelle fonction lesColonnes() et un exemple d'utilisation :
################################################################################################ lesColonnesV1 <- function(df,envar=FALSE,print=TRUE,order=TRUE) { ################################################################################################ nbc <- ncol(df) nbl <- nrow(df) cat("Voici les ",nbc,"colonnes (stat. sur ",nbl," lignes en tout)\n\n") noms <- names(df) nbvs <- rep(NA,nbc) mins <- rep(NA,nbc) maxs <- rep(NA,nbc) ddc <- cbind((1:nbc),nbvs,mins,maxs) row.names(ddc) <- noms colnames(ddc) <- c("Num","NbVal","Min","Max") pdc <- 1:nbc for (idc in pdc) { laCol <- as.numeric(df[,idc]) laCol <- laCol[!is.na(laCol)] ddc[idc,2] <- length(laCol) ddc[idc,3] <- min(laCol) ddc[idc,4] <- max(laCol) } # fin pour if (order) { # tri par ordre alphabétique des noms de variables idx <- order(noms) ddc <- ddc[idx,] } # fin si print.data.frame(as.data.frame(ddc),quote=FALSE,row.names=TRUE) cat("\n") if (envar) { return(as.data.frame(ddc)) } else { return() } # fin de si } # fin de fonction lesColonnesV1 ############################################################################## lesColonnesV2 <- function(df,envar=FALSE,print=TRUE,ordered=TRUE) { ################################################################# nbc <- ncol(df) nbl <- nrow(df) cat("Voici l'analyse des ",nbc,"colonnes (stat. sur ",nbl," lignes en tout)\n\n") noms <- names(df) nbvs <- rep(NA,nbc) mins <- rep(NA,nbc) maxs <- rep(NA,nbc) dist <- rep(NA,nbc) ddc <- cbind((1:nbc),nbvs,mins,maxs,dist) row.names(ddc) <- noms colnames(ddc) <- c("Num","NbVal","Min","Max","Distinct") pdc <- 1:nbc for (idc in pdc) { laCol <- as.numeric(df[,idc]) laCol <- laCol[!is.na(laCol)] ddc[idc,2] <- length(laCol) ddc[idc,3] <- min(laCol) ddc[idc,4] <- max(laCol) ddc[idc,5] <- length(unique(laCol)) } # fin pour if (ordered) { # tri par ordre alphabétique des noms de variables idx <- order(noms) ddc <- ddc[idx,] } # fin si print.data.frame(as.data.frame(ddc),quote=FALSE,row.names=TRUE) cat("\n") if (envar) { return(as.data.frame(ddc)) } else { return() } # fin de si } # fin de fonction lesColonnesV2 ################################################################# # démonstration : datagh("elf") cats("ancienne version") lesColonnesV1(elfdata) cats("nouvelle version") lesColonnesV2(elfdata)DOSSIER elf (après lecture de /home/gh/public_html/Datasets/elf.dar ) : Vous pouvez utiliser la matrice elfdata et les vecteurs ids(TXT) sexe (QLnum), sexe_(QLstr), age (QT), etud (QLnum), etud_ (QLstr) - les modalités sont nommées m_sexe m_etud - les unités sont nommées u_age - les intitulés sont nommées t_sexe t_age t_etud ancienne version ================ Voici les 6 colonnes (stat. sur 99 lignes en tout) Num NbVal Min Max AGE 2 99 11 78 ETUD 4 99 0 4 PROF 3 99 0 16 REGI 5 99 0 4 SEXE 1 99 0 1 USAG 6 99 0 3 NULL nouvelle version ================ Voici l'analyse des 6 colonnes (stat. sur 99 lignes en tout) Num NbVal Min Max valDistinct AGE 2 99 11 78 50 ETUD 4 99 0 4 5 PROF 3 99 0 16 13 REGI 5 99 0 4 5 SEXE 1 99 0 1 2 USAG 6 99 0 3 4 NULLCette fonction fait bien sur partie de statgh.r et se nomme lesColonnes().
Le problème est assez simple, puisqu'on veut seulement retirer des variables fortement corrélées. Deux boucles pour qui passent en revue la matrice triangulaire inférieure (ou supérieure) des coefficients de corrélation suffisent à repérer les variables à retirer, mises (via leur numéro de colonne) par notre programme dans le vecteur colasup. Une fois connues ces colonnes, on peut les retirer du jeu de données initial.
Il serait maladroit de vouloir imposer d'utiliser la corrélation de Kendall, ou de Spearman. Le mieux est de faire comme la fonction cor() : on laisse le choix de la méthode, avec la valeur "pearson" par défaut.
La notion de forte corrélation est très subjective. Utiliser le test de corrélation significative cor.test() est une solution un peu biaisée puisque ce test tend à trouver toute corrélation significative quand n est grand. Du coup, un seuil arbitraire mais variable et fixé par l'utilisateur est un choix raisonnable. Enfin, mettre les données réduites dans un nouveau fichier est sans doute une option utile, d'où les quatre paramètres utilisés par notre fonction...
################################################################# reduitVarsCor <- function(matdata,rho=0.8,file="",methode="pearson") { ################################################################# # on retire les variables fortement corrélées if (missing(matdata)) { cat(" reduitVars : élimine les colonnes linéairement corrélées à plus de rho, fourni en paramètre\n") cat(" (valeur par défaut : 0.8)\n") cat(" syntaxe : reduitVars(matdata,rho=0.8,file=\"\",methode=\"pearson\") \n\n") stop() } ; # fin si cat("\nRéduction du nombre de variables par suppression des variables fortement corrélées au seuil rho=",rho) cat("\npour la corrélation de ",methode,"\n\n") matdataorg <- matdata nbcol <- ncol(matdata) colasup <- c() # indices des colonnes à supprimer for (i in 1:(nbcol-1)) { vi <- matdata[,i] nomi <- sprintf("%-12s",names(matdata)[i]) for (j in (i+1):nbcol) { vj <- matdata[,j] corij <- cor(vi,vj,use="pairwise.complete.obs",method=methode) if (abs(corij)>=rho) { nomj <- sprintf("%-12s",names(matdata)[j]) cat(" la colonne ",sprintf("%5d",j)," soit ",nomj," est éliminée car fortement corrélée à ") cat("la colonne ",sprintf("%5d",i)," soit ",nomi," : rho = ",sprintf("%6.3f",corij),"\n") colasup <- c(colasup,j) } # fin de si } # fin pour j } # fin pour i newmatdata <- matdata cat(" réduction : ") nbcolasup <- length(colasup) if (nbcolasup>0) { nbcolasup <- length(colasup) if (nbcolasup==1) { cat(nbcolasup," colonne fortement corrélée linéairement supprimée \n") } else { cat(nbcolasup," colonnes fortement corrélées linéairement supprimées \n") } # fin si #print( colasup ) ###################################################### # # NE PAS ECRIRE # # colasup <- rev(sort(unique(colasup))) # # for (jcol in colasup) { # newmatdata <- newmatdata[,-jcol] # } # fin pour jcol # ###################################################### # # MAIS PLUTOT : # newmatdata <- newmatdata[,-colasup] } else { cat(" pas de colonnes fortement corrélées linéairement supprimées \n") } # fin si # recherche de la plus forte corrélation restante # affichage des dimensions avant et après cat("\nRésumé des dimensions\n") cat(" données initiales : ",sprintf("%5d",dim(matdataorg)),"\n") cat(" données finales : ",sprintf("%5d",dim(newmatdata)),"\n\n") nomi <- colnames(newmatdata)[1] nomj <- colnames(newmatdata)[1] maxcor <- 0 nbcol <- ncol(newmatdata) for (i in 1:(nbcol-1)) { vi <- newmatdata[,i] for (j in (i+1):nbcol) { vj <- newmatdata[,j] corij <- cor(vi,vj,use="pairwise.complete.obs",method=methode) if (abs(corij)>=abs(maxcor)) { maxcor <- corij imax <- i jmax <- j nomi <- colnames(newmatdata)[i] nomj <- colnames(newmatdata)[j] } # fin de si } # fin pour j } # fin pour i cat(" le plus fort coefficient de corrélation restant est ",sprintf("%6.3f",maxcor),"\n") cat(" et met en jeu les colonnes ",imax," et ",jmax," soit ",nomi," et ",nomj,"\n\n") if (file!="") { write.table(newmatdata,file=file,quote=FALSE,) cat("les données finales ont été écrites dans le fichier ",file,"\n") } ; # fin si return(newmatdata) } # fin de fonction reduitVarsCor ################################################################# # exemple d'utilisation : lea <- lit.dar("leadb710x46.dar.txt")[ ,-1] # sans la colonne 1 ! leaReduit <- reduitVarsCor(matdata=lea,rho=0.8,methode=methode)Réduction du nombre de variables par suppression des variables fortement corrélées au seuil rho= 0.8 pour la corrélation de pearson la colonne 3 soit MW est éliminée car fortement corrélée à la colonne 1 soit Length : rho = 0.995 la colonne 6 soit Netcharge est éliminée car fortement corrélée à la colonne 2 soit pI : rho = 0.923 la colonne 38 soit D.E.K.R.1 est éliminée car fortement corrélée à la colonne 2 soit pI : rho = -0.911 la colonne 5 soit Gravy est éliminée car fortement corrélée à la colonne 4 soit Foldindex : rho = 0.987 la colonne 8 soit Hydrophobi est éliminée car fortement corrélée à la colonne 4 soit Foldindex : rho = 0.879 la colonne 11 soit Buried est éliminée car fortement corrélée à la colonne 4 soit Foldindex : rho = 0.811 la colonne 13 soit Transmembr est éliminée car fortement corrélée à la colonne 4 soit Foldindex : rho = 0.905 la colonne 8 soit Hydrophobi est éliminée car fortement corrélée à la colonne 5 soit Gravy : rho = 0.873 la colonne 13 soit Transmembr est éliminée car fortement corrélée à la colonne 5 soit Gravy : rho = 0.887 la colonne 38 soit D.E.K.R.1 est éliminée car fortement corrélée à la colonne 6 soit Netcharge : rho = -0.997 la colonne 8 soit Hydrophobi est éliminée car fortement corrélée à la colonne 7 soit Hydrophili : rho = -0.896 la colonne 13 soit Transmembr est éliminée car fortement corrélée à la colonne 7 soit Hydrophili : rho = -0.936 la colonne 17 soit E.Unip est éliminée car fortement corrélée à la colonne 7 soit Hydrophili : rho = 0.855 la colonne 35 soit D.E est éliminée car fortement corrélée à la colonne 7 soit Hydrophili : rho = 0.856 la colonne 36 soit K.R est éliminée car fortement corrélée à la colonne 7 soit Hydrophili : rho = 0.850 la colonne 37 soit D.E.K.R est éliminée car fortement corrélée à la colonne 7 soit Hydrophili : rho = 0.918 la colonne 11 soit Buried est éliminée car fortement corrélée à la colonne 8 soit Hydrophobi : rho = 0.832 la colonne 13 soit Transmembr est éliminée car fortement corrélée à la colonne 8 soit Hydrophobi : rho = 0.926 la colonne 39 soit A.I.L.V est éliminée car fortement corrélée à la colonne 9 soit Flexibility : rho = -0.828 la colonne 15 soit G.Unip est éliminée car fortement corrélée à la colonne 10 soit Bulkiness : rho = -0.882 la colonne 13 soit Transmembr est éliminée car fortement corrélée à la colonne 11 soit Buried : rho = 0.868 la colonne 35 soit D.E est éliminée car fortement corrélée à la colonne 17 soit E.Unip : rho = 0.864 la colonne 37 soit D.E.K.R est éliminée car fortement corrélée à la colonne 17 soit E.Unip : rho = 0.870 la colonne 37 soit D.E.K.R est éliminée car fortement corrélée à la colonne 19 soit K.Unip : rho = 0.810 la colonne 42 soit S.T est éliminée car fortement corrélée à la colonne 25 soit T.Unip : rho = 0.807 la colonne 41 soit N.Q est éliminée car fortement corrélée à la colonne 32 soit Q.Unip : rho = 0.903 la colonne 37 soit D.E.K.R est éliminée car fortement corrélée à la colonne 35 soit D.E : rho = 0.949 la colonne 37 soit D.E.K.R est éliminée car fortement corrélée à la colonne 36 soit K.R : rho = 0.904 la colonne 45 soit C.F.Y.W est éliminée car fortement corrélée à la colonne 40 soit F.W.Y : rho = 0.967 réduction : 16 colonnes fortement corrélées linéairement supprimées Résumé des dimensions données initiales : 710 45 données finales : 710 29 le plus fort coefficient de corrélation restant est 0.796 et met en jeu les colonnes 4 et 11 soit Hydrophili et K.UnipLe problème n'est pas simple, et en tous cas beaucoup plus compliqué qu'à la question précédente. En effet, si X et Y sont corrélées, il est possible que Y soit aussi corrélée à Z sans que X ne soit corrélée à Z car il n'y a pas transitivité, comme on dit en mathématique, de la notion de corrélation. Une solution consiste, à partir d'une variable, lui adjoindre toutes celles qui lui sont fortement corrélées, puis à calculer la clotûre transitive, c'est-à-dire ajouter aussi les autres variables qui sont fortement corrélées aux variables qu'on vient d'ajouter et de recommencer tant qu'on ajoute des variables...
Dans le programme qui suit, on utilise 2 vecteurs nommés colatrait et colvues. Au départ, colatrait vaut 1 pour toutes les variables et on met 0 dans colatrait[j] quand on traité la variable j. Il y a donc une boucle principale basée sur le nombre de variables encore à 1 dans colatrait.
Le tableau colvues permet de savoir quel est l'état de la variable : déjà traitée, en cours ou non encore vue. Afin de savoir quelles variables sont directement corrélées et celles qui sont ajoutées ensuite, notre fonction insère une étoile pour la corrélation directe et deux étoiles pour la corrélation transitive.
################################################################################### corClust <- function(data,seuil=0.8,methode="pearson",detail=FALSE) { ################################################################################### # on regroupe transitivement les variables fortement corrélées if (missing(data)) { cat(" corClust : regroupe les colonnes linéairement corrélées pour rho supérieur au seuil fourni en paramètre\n") cat(" (valeur par défaut : 0.8)\n") cat(" syntaxe : corClust(matdata,seuil=0.8,methode=\"pearson\") \n\n") stop() } ; # fin si matdata <- data rho <- seuil cat("Regroupement transitif de variables fortement corrélées au seuil rho=",rho,"\n") cat("pour la corrélation de ",methode,"\n") nbcol <- ncol(matdata) nomcol <- names(matdata) colatrait <- rep(1,nbcol) colvues <- rep(2,nbcol) # 2 : colonne non vue, 1 en cours, 0 déja traitée lesgrp <- rep(0,nbcol) # numéro du groupe # on passe les colonnes en revue : celles qui sont à rho>seuil sont mises # ensemble (une étoile) ; on met alors colatrait (colonnes à traiter) à 0 et on # recommence avec les variables du groupe par transitivité (deux étoiles) jusqu'à saturation # puis on passe au groupe suivant ... nb1 <- sum(colatrait==1) grp <- 0 # numéro de groupe while (nb1>0) { ibase <- 1 while (ibase<=nbcol) { if (colatrait[ibase]==1) { itrait <- ibase ibase <- nbcol + 1 } # fin de si sur colonne à traiter ibase <- ibase + 1 } # fin tant que sur ibase grp <- grp + 1 eff <- 1 vi <- matdata[,itrait] colvues[itrait] <- 0 cat("\nGroupe ",sprintf("%2d",grp)," : ") cat(" ",nomcol[itrait]) lesgrp[itrait] <- grp if (itrait<=nbcol) { colatrait[itrait] <- 0 # premier passage eff <- 0 for (j in (itrait+1):nbcol) { vj <- matdata[,j] corij <- cor(vi,vj,use="pairwise.complete.obs",method=methode) if (abs(corij)>=rho) { colvues[j] <- 1 cat(" *",nomcol[j]) lesgrp[j] <- grp eff <- eff + 1 } # fin de si } # fin pour j # deuxieme passage (transitivité faible de rho) nb2 <- sum(colvues==1) jtrait <- nbcol+1 while (nb2>0) { jbase <- 1 while (jbase<=nbcol) { if ( (colvues[jbase]==1) & (colatrait[jbase]==1) ) { jtrait <- jbase jbase <- nbcol + 1 colvues[jtrait] <- 0 colatrait[jtrait] <- 0 } # fin de si sur colonne à traiter jbase <- jbase + 1 } # fin tant que sur jbase nb3 <- sum(colvues==1) if (jtrait<=nbcol) { vjt <- matdata[,jtrait] for (j in (1:nbcol)) { if (colvues[j]==2) { if (colatrait[j]==1) { vj <- matdata[,j] corij <- cor(vjt,vj,use="pairwise.complete.obs",method=methode) if (abs(corij)>=rho) { colvues[j] <- 1 cat(" **",nomcol[j]) lesgrp[j] <- grp } # fin de si } # fin de si } # fin de si } # fin pour j } # fin si nb2 <- sum(colvues==1) } # fin de tant que } # fin de si on a au moins une colonne à traiter nb1 <- sum(colatrait==1) } # fin de tant que cat("\n") # mise en forme des résultats retenue <- rep(0,nbcol) gard <- rep(0,grp) for (igrp in (1:grp)) { icol <- 1 igard <- 0 while (icol<=nbcol) { if (lesgrp[icol]==igrp) { igard <- icol icol <- nbcol + 1 } # fin si icol <- icol + 1 } # fin tant que gard[igrp] <- igard retenue[igard] <- igrp } # fin pour igrp mret <- cbind(lesgrp,retenue) row.names(mret) <- names(matdata) colnames(mret) <- c("Groupe","Variable retenue") if (detail) { cat("\n") cat("Numéros de groupe de corrélation et variable retenue : \n") print(mret) cat("\n") cat("On retient ",grp," variables :") } # fin si garde <- matdata[,gard] cat("\nRésumé des dimensions\n") cat(" données initiales : ",sprintf("%5d",dim(matdata)),"\n") cat(" données finales : ",sprintf("%5d",dim(garde)),"\n\n") if (detail) { return(garde) } # fin si } # fin de fonction corClust ################################################################# # exemple d'utilisation : lea <- lit.dar("leadb710x46.dar.txt")[ ,-1] # sans la colonne 1 ! corClust(data=lea,seuil=0.8,methode="pearson",detail=TRUE)Regroupement transitif de variables fortement corrélées au seuil rho= 0.8 pour la corrélation de pearson Groupe 1 : Length * MW Groupe 2 : pI * Netcharge * D.E.K.R.1 Groupe 3 : Foldindex * Gravy * Hydrophobi * Buried * Transmembr ** Hydrophili ** E.Unip ** D.E ** K.R ** D.E.K.R ** K.Unip Groupe 4 : Flexibility * A.I.L.V Groupe 5 : Bulkiness * G.Unip Groupe 6 : Accres Groupe 7 : MW.Length Groupe 8 : D.Unip Groupe 9 : R.Unip Groupe 10 : A.Unip Groupe 11 : L.Unip Groupe 12 : I.Unip Groupe 13 : V.Unip Groupe 14 : S.Unip Groupe 15 : T.Unip * S.T Groupe 16 : C.Unip Groupe 17 : M.Unip Groupe 18 : F.Unip Groupe 19 : W.Unip Groupe 20 : Y.Unip Groupe 21 : N.Unip Groupe 22 : Q.Unip * N.Q Groupe 23 : P.Unip Groupe 24 : H.Unip Groupe 25 : F.W.Y * C.F.Y.W Groupe 26 : C.W Groupe 27 : R.E.S.P Numéros de groupe de corrélation et variable retenue : Groupe Variable retenue Length 1 1 pI 2 2 MW 1 0 Foldindex 3 3 Gravy 3 0 Netcharge 2 0 Hydrophili 3 0 Hydrophobi 3 0 Flexibility 4 4 Bulkiness 5 5 Buried 3 0 Accres 6 6 Transmembr 3 0 MW.Length 7 7 G.Unip 5 0 D.Unip 8 8 E.Unip 3 0 R.Unip 9 9 K.Unip 3 0 A.Unip 10 10 L.Unip 11 11 I.Unip 12 12 V.Unip 13 13 S.Unip 14 14 T.Unip 15 15 C.Unip 16 16 M.Unip 17 17 F.Unip 18 18 W.Unip 19 19 Y.Unip 20 20 N.Unip 21 21 Q.Unip 22 22 P.Unip 23 23 H.Unip 24 24 D.E 3 0 K.R 3 0 D.E.K.R 3 0 D.E.K.R.1 2 0 A.I.L.V 4 0 F.W.Y 25 25 N.Q 22 0 S.T 15 0 C.W 26 26 R.E.S.P 27 27 C.F.Y.W 25 0 On retient 27 variables : Résumé des dimensions données initiales : 710 45 données finales : 710 27 [... données non affichées]Pour comprendre ce qu'a fait le programme, on peut calculer la matrice des corrélations des variables du groupe 3 :
lea <- lit.dar("leadb710x46.dar.txt")[ ,-1] # sans la colonne 1 ! # démonstration avec le groupe 3 nomvar <- c("Foldindex","Gravy","Hydrophobi","Buried","Transmembr","Hydrophili","E.Unip","D.E","K.R","D.E.K.R","K.Unip") varcor <- lea[ , nomvar ] mdc(varcor)Matrice des corrélations au sens de Pearson pour 710 lignes et 11 colonnes Foldindex Gravy Hydrophobi Buried Transmembr Hydrophili E.Unip D.E K.R D.E.K.R K.Unip Foldindex 1.000 Gravy 0.987 1.000 Hydrophobi 0.879 0.873 1.000 Buried 0.811 0.774 0.832 1.000 Transmembr 0.905 0.887 0.926 0.868 1.000 Hydrophili -0.748 -0.719 -0.896 -0.766 -0.936 1.000 E.Unip -0.541 -0.490 -0.678 -0.664 -0.738 0.855 1.000 D.E -0.446 -0.394 -0.605 -0.633 -0.721 0.856 0.864 1.000 K.R -0.477 -0.439 -0.768 -0.670 -0.716 0.850 0.735 0.723 1.000 D.E.K.R -0.494 -0.444 -0.725 -0.698 -0.773 0.918 0.870 0.949 0.904 1.000 K.Unip -0.536 -0.494 -0.628 -0.729 -0.786 0.796 0.682 0.723 0.793 0.810 1.000La variable à traiter au départ était Foldindex. Le premier passage est venu ajouter les 4 variables Gravy Hydrophobi Buried Transmembr car leur rho avec Foldindex est supérieur à 0,8. Ensuite, la variable à traiter est Gravy mais elle n'ajoute aucune variable. Puis on traite Hydrophobi qui ajoute Hydrophili...
Il est possible d'ajouter une variable en paramètre pour afficher le suivi des ajouts de corrélation :
################################################################################### corClust <- function(data,seuil=0.8,methode="pearson",detail=FALSE,suivi=FALSE) { ################################################################################### # on regroupe transitivement les variables fortement corrélées if (missing(data)) { cat(" corClust : regroupe les colonnes linéairement corrélées pour rho supérieur au seuil fourni en paramètre\n") cat(" (valeur par défaut : 0.8)\n") cat(" syntaxe : corClust(matdata,seuil=0.8,methode=\"pearson\") \n\n") stop() } ; # fin si matdata <- data rho <- seuil cat("Regroupement transitif de variables fortement corrélées au seuil rho=",rho,"\n") cat("pour la corrélation de ",methode,"\n") nbcol <- ncol(matdata) nomcol <- names(matdata) colatrait <- rep(1,nbcol) colvues <- rep(2,nbcol) # 2 : colonne non vue, 1 en cours, 0 déja traitée lesgrp <- rep(0,nbcol) # numéro du groupe # on passe les colonnes en revue : celles qui sont à rho>seuil sont mises # ensemble (une étoile) ; on met alors colatrait (colonnes à traiter) à 0 et on # recommence avec les variables du groupe par transitivité (deux étoiles) jusqu'à saturation # puis on passe au groupe suivant ... nb1 <- sum(colatrait==1) grp <- 0 # numéro de groupe while (nb1>0) { ibase <- 1 while (ibase<=nbcol) { if (colatrait[ibase]==1) { itrait <- ibase ibase <- nbcol + 1 } # fin de si sur colonne à traiter ibase <- ibase + 1 } # fin tant que sur ibase grp <- grp + 1 eff <- 1 vi <- matdata[,itrait] colvues[itrait] <- 0 cat("\nGroupe ",sprintf("%2d",grp)," : ") if (suivi) { cat(" basé sur ") } cat(" ",nomcol[itrait]) if (suivi) { cat("\n") } lesgrp[itrait] <- grp if (itrait<=nbcol) { colatrait[itrait] <- 0 # premier passage if (suivi) { cat(" variables corrélées directement : ") } eff <- 0 for (j in (itrait+1):nbcol) { vj <- matdata[,j] corij <- cor(vi,vj,use="pairwise.complete.obs",method=methode) if (abs(corij)>=rho) { colvues[j] <- 1 cat(" *",nomcol[j]) lesgrp[j] <- grp eff <- eff + 1 } # fin de si } # fin pour j if (suivi) { cat("\n") } # deuxieme passage (transitivité faible de rho) if (suivi) { cat(" variables corrélées transitivement : \n") } nb2 <- sum(colvues==1) jtrait <- nbcol+1 while (nb2>0) { jbase <- 1 while (jbase<=nbcol) { if ( (colvues[jbase]==1) & (colatrait[jbase]==1) ) { jtrait <- jbase jbase <- nbcol + 1 colvues[jtrait] <- 0 colatrait[jtrait] <- 0 } # fin de si sur colonne à traiter jbase <- jbase + 1 } # fin tant que sur jbase nb3 <- sum(colvues==1) if (jtrait<=nbcol) { vjt <- matdata[,jtrait] for (j in (1:nbcol)) { if (colvues[j]==2) { if (colatrait[j]==1) { vj <- matdata[,j] corij <- cor(vjt,vj,use="pairwise.complete.obs",method=methode) if (abs(corij)>=rho) { colvues[j] <- 1 if (suivi) { cat(" ",nomcol[jtrait]) } cat(" **",nomcol[j]) if (suivi) { cat("\n") } lesgrp[j] <- grp } # fin de si } # fin de si } # fin de si } # fin pour j } # fin si nb2 <- sum(colvues==1) } # fin de tant que if (suivi) { cat("\n") } } # fin de si on a au moins une colonne à traiter nb1 <- sum(colatrait==1) # pour debug : attends() } # fin de tant que cat("\n") # mise en forme des résultats retenue <- rep(0,nbcol) gard <- rep(0,grp) for (igrp in (1:grp)) { icol <- 1 igard <- 0 while (icol<=nbcol) { if (lesgrp[icol]==igrp) { igard <- icol icol <- nbcol + 1 } # fin si icol <- icol + 1 } # fin tant que gard[igrp] <- igard retenue[igard] <- igrp } # fin pour igrp mret <- cbind(lesgrp,retenue) row.names(mret) <- names(matdata) colnames(mret) <- c("Groupe","Variable retenue") if (detail) { cat("\n") cat("Numéros de groupe de corrélation et variable retenue : \n") print(mret) cat("\n") cat("On retient ",grp," variables :") } # fin si garde <- matdata[,gard] cat("\nRésumé des dimensions\n") cat(" données initiales : ",sprintf("%5d",dim(matdata)),"\n") cat(" données finales : ",sprintf("%5d",dim(garde)),"\n\n") if (detail) { return(garde) } # fin si } # fin de fonction corClust ################################################################# # exemple d'utilisation : lea <- lit.dar("leadb710x46.dar.txt")[ ,-1] # sans la colonne 1 ! leaReduit <- corClust(data=lea,seuil=0.8,methode="pearson",suivi=TRUE)Regroupement transitif de variables fortement corrélées au seuil rho= 0.8 pour la corrélation de pearson Groupe 1 : basé sur Length variables corrélées directement : * MW variables corrélées transitivement : Groupe 2 : basé sur pI variables corrélées directement : * Netcharge * D.E.K.R.1 variables corrélées transitivement : Groupe 3 : basé sur Foldindex variables corrélées directement : * Gravy * Hydrophobi * Buried * Transmembr variables corrélées transitivement : Hydrophobi ** Hydrophili Hydrophili ** E.Unip Hydrophili ** D.E Hydrophili ** K.R Hydrophili ** D.E.K.R D.E.K.R ** K.Unip Groupe 4 : basé sur Flexibility variables corrélées directement : * A.I.L.V variables corrélées transitivement : Groupe 5 : basé sur Bulkiness variables corrélées directement : * G.Unip variables corrélées transitivement : Groupe 6 : basé sur Accres variables corrélées directement : variables corrélées transitivement : Groupe 7 : basé sur MW.Length variables corrélées directement : variables corrélées transitivement : Groupe 8 : basé sur D.Unip variables corrélées directement : variables corrélées transitivement : Groupe 9 : basé sur R.Unip variables corrélées directement : variables corrélées transitivement : Groupe 10 : basé sur A.Unip variables corrélées directement : variables corrélées transitivement : Groupe 11 : basé sur L.Unip variables corrélées directement : variables corrélées transitivement : Groupe 12 : basé sur I.Unip variables corrélées directement : variables corrélées transitivement : Groupe 13 : basé sur V.Unip variables corrélées directement : variables corrélées transitivement : Groupe 14 : basé sur S.Unip variables corrélées directement : variables corrélées transitivement : Groupe 15 : basé sur T.Unip variables corrélées directement : * S.T variables corrélées transitivement : Groupe 16 : basé sur C.Unip variables corrélées directement : variables corrélées transitivement : Groupe 17 : basé sur M.Unip variables corrélées directement : variables corrélées transitivement : Groupe 18 : basé sur F.Unip variables corrélées directement : variables corrélées transitivement : Groupe 19 : basé sur W.Unip variables corrélées directement : variables corrélées transitivement : Groupe 20 : basé sur Y.Unip variables corrélées directement : variables corrélées transitivement : Groupe 21 : basé sur N.Unip variables corrélées directement : variables corrélées transitivement : Groupe 22 : basé sur Q.Unip variables corrélées directement : * N.Q variables corrélées transitivement : Groupe 23 : basé sur P.Unip variables corrélées directement : variables corrélées transitivement : Groupe 24 : basé sur H.Unip variables corrélées directement : variables corrélées transitivement : Groupe 25 : basé sur F.W.Y variables corrélées directement : * C.F.Y.W variables corrélées transitivement : Groupe 26 : basé sur C.W variables corrélées directement : variables corrélées transitivement : Groupe 27 : basé sur R.E.S.P variables corrélées directement : variables corrélées transitivement : Résumé des dimensions données initiales : 710 45 données finales : 710 27On remarquera que cette technique est plus efficace que la précédente pour l'exemple considéré car elle ne conserve que 27 variables au lieu de 29.
Trouver le meilleur représentant par paquet n'est pas un problème simple non plus car il faut savoir qu'une variable "fortement explicative" lorsqu'on analyse tous les modèles à une seule variable ne sera pas forcément sélectionnée lorsqu'on passe en revue tous les modèles à deux, trois variables ou plus...
La réponse est Rweb. Il s'agit d'une interface Web à R sous la forme d'un formulaire. En France, on peut par exemple utiliser pbil-rweb.
En général, Rweb contient de nombreux packages. Vous pouvez par exemple essayer de copier/coller le code suivant pour le vérifier :
## récupération d'un tableau XHTML dans R # chargement du package XML library(XML) # chargement du package stringr library(stringr) # pour la fonction str_replace() # url de la page Web à utiliser url <- "http://forge.info.univ-angers.fr/~gh/wstat/Programmation_R/progr.php?n=1&m=s" tabp <- readHTMLTable(url,which=3,stringsAsFactors=FALSE) print(tabp,right=FALSE) # un peu de "ménage" colnames(tabp) <- str_replace(colnames(tabp),"\\n","") # tout sauf les deux derniers caractères tsdc <- function(x) { substr(x,1,nchar(x)-2) } tabp[,2] <- as.numeric( tsdc(tabp[,2]) ) print(tabp,right=FALSE) ## récupération d'une liste ip <- "http://fr.wikipedia.org/wiki/Institut_Pasteur" cdr <- readHTMLList(ip,which=7) # conversion des accents iconv(cdr,from="UTF-8",to="LATIN1") print(cbind(cdr))Oui, Sweave est indispensable pour les rapports et les articles car c'est un gage de reproductibilité et de gain de temps. Après tout, il ne s'agit que d'inclure du R dans du code LaTeX. Ce n'est donc pas compliqué quand on connait LaTeX. Voici un extrait de code Sweave :
\section{Introduction} <<echo=FALSE,results=hide>>= options(encoding = "latin1") source("poursweave.r",encoding="latin1") nlu <- 10 # nombre de lignes \ufffd lire usr <- "INRA" # pour montrer que c'est interactif dos <- "iris" # iris ou elf nbl <- demoSw(dos,nlu,usr,1) @ \section{Quelques lignes de donn\'{e}es du fichier : \Sexpr{dos} } Le fichier elf contient \Sexpr{nbl} lignes. L'utilisateur \Sexpr{usr} veut traiter les \Sexpr{nlu} premi\`{e}res lignes de ce fichier.Comme on peut le voir, les fragments \Sexpr{...} correspondent à des valeurs que R insère dans le code LaTeX alors que les parties qui commencent par << et qui se terminent par @ sont du code R directement interprété.
Oui, interfacer R et php est assez facile à réaliser car R a un mode batch. En d'autres termes, on peut au moins communiquer via des fichiers-texte.
En voici un exemple, nommé aqt.