Programmer en R Session de formation, décembre 2013 gilles.hunault "at" univ-angers.fr                        Solutions pour la séance numéro 4 (énoncés) Avec debug(f), on dispose d'une exécution en pas à pas de la fonction f ; les commandes sont n (ou la touche Entrée) pour next, c'est-à-dire poursuivre l'exécution de l'instruction, c pour continuer, c'est-à-dire passer à l'instruction suivante et Q pour quitter. A l'intérieur d'une boucle, n passe à l'itération suivante alors que c exécute tout la boucle. On dispose aussi de where qui indique d'où provient l'exécution (lignes de fichiers-sources, fonctions appelantes...). Pour mettre fin au mode de déboggage de la fonction f on peut taper undebug(f) (en dehors de la session de débogage) ou recharger la définition de la fonction avec source(...) si la fonction est définie dans un fichier. En mode debug, la fonction browser() est exécutée, qui permet l'inspection de l'environnement. On peut l'appeler directement quand on veut avoir le détail des valeurs des variables et de l'environnement. Taper n, c ou Q vient exécuter les commandes indiquées précédemment. Si on a nommé une variable n, il faut écrire explicitement print(n) pour connaitre la valeur de n. En mode browser, on peut exécuter toute commande R pour sauvegarder des variables, pour tester une expression... Voici un exemple de session pour le cas de la fonction sdu() : > sdu( c(1,0,0,1,1,1,0,1,1),2 ) # fonction incorrecte !      [1] 3 4 5 7            > debug(sdu)            #######################################################            > sdu( c(1,0,0,1,1,1,0,1,1),2 ) # on relance l'exécution (en mode debug)            debugging in: sdu(c(1, 0, 0, 1, 1, 1, 0, 1, 1), 2)      debug à sdu.r#5 :{      n <- length(x)      runs <- NULL      for (i in 1:(n - k)) {      if (all(x[i:i + k - 1] == 1)) {      runs <- c(runs, i)      }      }      return(runs)      }            Browse[2]> n # taper n fait "next"            debug à sdu.r#7 :n <- length(x)            Browse[2]> n            debug à sdu.r#8 :runs <- NULL            Browse[2]> i # l'objet i n'est pas encore défini            Erreur : objet 'i' introuvable            Browse[2]> n            debug à sdu.r#9 :for (i in 1:(n - k)) {      if (all(x[i:i + k - 1] == 1)) {      runs <- c(runs, i)      }      }            Browse[2]> n            debug à sdu.r#9 :i            Browse[2]> i            NULL            Browse[2]> n            debug à sdu.r#10 :if (all(x[i:i + k - 1] == 1)) {      runs <- c(runs, i)      }            Browse[2]> n            debug à sdu.r#10 :NULL            Browse[2]> n            debug à sdu.r#9 :i            Browse[2]> i            [1] 1            Browse[2]> n            debug à sdu.r#10 :if (all(x[i:i + k - 1] == 1)) {      runs <- c(runs, i)      }            Browse[2]> n            debug à sdu.r#10 :NULL            Browse[2]> n            debug à sdu.r#9 :i            Browse[2]> n            debug à sdu.r#10 :if (all(x[i:i + k - 1] == 1)) {      runs <- c(runs, i)      }            Browse[2]> i      [1] 3            Browse[2]> i+k-1      [1] 4            Browse[2]> i:i+k-1      [1] 4            Browse[2]> i:(i+k-1)      [1] 3 4            Browse[2]> n      debug à sdu.r#11 :runs <- c(runs, i)            Browse[2]> i      [1] 3            Browse[2]> c      debug à sdu.r#14 :return(runs)            Browse[2]> i      [1] 7            Browse[2]> n      exiting from: sdu(c(1, 0, 0, 1, 1, 1, 0, 1, 1), 2)      [1] 3 4 5 7            > i      Erreur : objet 'i' introuvable            #######################################################            > undebug(sdu)                  On en déduit que l'erreur vient de la plage de valeurs x[ i : (i+k-1) ] incorrectement écrite en x[ i : i+k-1 ]. Voici donc le bon code pour la fonction findruns() : findruns <- function(x,k) {            n <- length(x)      runs <- NULL      for (i in 1:(n-k+1)) {      if (all(x[i:(i+k-1)]==1)) {      runs <- c(runs,i)      } # fin si      } # fin pour i      return(runs)      } # fin de fonction findruns            print( findruns( c(1,0,0,1,1,1,0,1,1),2 )) # affiche bien 4 5 8                  Rappel : une fois la séance de debug terminée, il faut écrire undebug(f) pour exécuter normalement la fonction f ou recharger sa définition. Pour avoir un retour sur l'erreur, le mieux est d'utiliser options(error=recover) puis de repasser à options(error=NULL) quand tout est réglé. Démonstration : > sinksrc("cites.r") # en mode normal, une erreur quitte l'évaluation du code            Erreur dans d[i, j] : indice hors limites            #######################################################            > options(error=utils::recover) # désormais on pourra revenir à l'exécution            #######################################################            > sinksrc("cites.r")            Erreur dans mind(matq) : indice hors limites            Enter a frame number, or 0 to exit            1: sinksrc("cites.r")      2: statgh.r#9243: sinksource(fn, encoding = "latin1")      3: statgh.r#9211: source(fp, encoding = encoding, print.eval = TRUE)      4: statgh.r#9211: eval.with.vis(ei, envir)      5: statgh.r#9211: eval.with.vis(expr, envir, enclos)      6: statgh.r#9211: mind(matq)            Sélection : 6      Called from: eval(substitute(browser(skipCalls = skip), list(skip = 7 - which)),      envir = sys.frame(which))            Browse[1]> d      [,1] [,2] [,3] [,4] [,5]      [1,] 0 12 13 8 20      [2,] 12 0 15 28 88      [3,] 13 15 0 6 9      [4,] 8 28 6 0 33      [5,] 20 88 9 33 0            Browse[1]> nb      [1] 5            Browse[1]> dd      [,1] [,2] [,3] [,4] [,5] [,6]      [1,] 0 12 13 8 20 1      [2,] 12 0 15 28 88 2      [3,] 13 15 0 6 9 3      [4,] 8 28 6 0 33 4      [5,] 20 88 9 33 0 5            Browse[1]> wmins      [,1] [,2] [,3] [,4]      [1,] 4 3 4 5      [2,] 8 15 6 33            Browse[1]> i      [1] 2            Browse[1]> j      [1] 15            Browse[1]> d[i,j]            Erreur pendant l'emballage (wrapup) : indice hors limites      Browse[1]>Q            #######################################################            > options(error=NULL) # il est prudent de remettre l'option à sa valeur par défaut                  On peut voir qu'il y a deux erreurs, basées sur la confusion entre les dimensions 1 (les lignes) et 2 (les colonnes). Il faut donc écrire i <- which.min(wmins[2,]) au lieu de i <- which.min(wmins[1,]) De même, il faut remplacer j <- wmins[1,i] par j <- wmins[2,i] Voici donc le bon code pour la fonction mind(), qui renvoie bien 6 3 4 ce qui doit s'interpréter comme distance minimale 6 vue en colonne 3, ligne 4. mind <- function( d ) {            imin <- function(x) { # cette fonction est "locale" à mind      lx <- length(x)      i <- x[lx]      j <- which.min(x[(i+1):(lx-1)])      k <- i+j            return( c(k,x[k]) )            } # fin de fonction imin            nb <- nrow(d)      dd <- cbind(d,1:nb) # rappel du numéro de ligne      wmins <- apply(dd[-nb,],1,imin)      i <- which.min(wmins[2,])      j <- wmins[1,i]            return( c(d[i,j],i,j) )            } # fin de fonction mind            matq <- matrix(nrow=5,ncol=5)      matq[1,] <- c( 0,12,13, 8,20)      matq[2,] <- c(12, 0,15,28,88)      matq[3,] <- c(13,15, 0, 6, 9)      matq[4,] <- c( 8,28, 6, 0,33)      matq[5,] <- c(20,88, 9,33, 0)            dm <- mind( matq )      cat(" distance min :",dm[1]," vue en col lig ",dm[-1],"\n")                  Les fonctions pour le profilage sont Rprof() et summaryRprof(). Voici ce que cela donne sur notre exemple. Il faut aussi consulter le fichier Rprof.out généré par R. Rprof()      sinksrc("pdm.r")      Rprof(NULL)      summaryRprof()                  powers1      =======            prod prod prod prod      [1,] 1 1 1 1 1      [2,] 2 4 8 16 32      [3,] 3 9 27 81 243      [4,] 4 16 64 256 1024      [5,] 5 25 125 625 3125            powers2      =======            [,1] [,2] [,3] [,4] [,5]      [1,] 1 1 1 1 1      [2,] 2 4 8 16 32      [3,] 3 9 27 81 243      [4,] 4 16 64 256 1024      [5,] 5 25 125 625 3125            les 4 méthodes      ==============            prod prod prod prod      [1,] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1      [2,] 2 4 8 16 32 2 4 8 16 32 2 4 8 16 32 2 4 8 16 32      [3,] 3 9 27 81 243 3 9 27 81 243 3 9 27 81 243 3 9 27 81 243      [4,] 4 16 64 256 1024 4 16 64 256 1024 4 16 64 256 1024 4 16 64 256 1024      [5,] 5 25 125 625 3125 5 25 125 625 3125 5 25 125 625 3125 5 25 125 625 3125            comparaison des durées des 4 méthodes      =====================================            powers1 durée = 15.34975 s      powers2 durée = 0.1522987 s      powers3 durée = 0.8030689 s      powers4 durée = 0.763536 s                  user.self sys.self elapsed user.child sys.child      t1 12.41 2.85 15.281 0 0      t2 0.09 0.01 0.100 0 0      t3 0.70 0.04 0.749 0 0      t4 0.65 0.04 0.709 0 0                  vous pouvez utiliser pdm.sor                  $by.self      self.time self.pct total.time total.pct      "cbind" 15.20 88.89 15.20 88.89      "outer" 0.76 4.44 0.76 4.44      "apply" 0.34 1.99 0.54 3.16      "gc" 0.24 1.40 0.24 1.40      "matrix" 0.10 0.58 0.10 0.58      "*" 0.08 0.47 0.08 0.47      "aperm.default" 0.08 0.47 0.08 0.47      "t.default" 0.08 0.47 0.08 0.47      "powers2" 0.06 0.35 0.10 0.58      "array" 0.06 0.35 0.06 0.35      "unlist" 0.06 0.35 0.06 0.35      "powers1" 0.02 0.12 15.28 89.36      "print.default" 0.02 0.12 0.02 0.12            $by.total      total.time total.pct self.time self.pct      "eval.with.vis" 17.10 100.00 0.00 0.00      "sinksource" 17.10 100.00 0.00 0.00      "sinksrc" 17.10 100.00 0.00 0.00      "source" 17.10 100.00 0.00 0.00      "duree" 17.08 99.88 0.00 0.00      "eval" 17.08 99.88 0.00 0.00      "eval.parent" 17.08 99.88 0.00 0.00      "system.time" 17.08 99.88 0.00 0.00      "powers1" 15.28 89.36 0.02 0.12      "cbind" 15.20 88.89 15.20 88.89      "outer" 0.76 4.44 0.76 4.44      "powers3" 0.76 4.44 0.00 0.00      "powers4" 0.70 4.09 0.00 0.00      "t" 0.62 3.63 0.00 0.00      "apply" 0.54 3.16 0.34 1.99      "gc" 0.24 1.40 0.24 1.40      "matrix" 0.10 0.58 0.10 0.58      "powers2" 0.10 0.58 0.06 0.35      "*" 0.08 0.47 0.08 0.47      "aperm.default" 0.08 0.47 0.08 0.47      "t.default" 0.08 0.47 0.08 0.47      "aperm" 0.08 0.47 0.00 0.00      "array" 0.06 0.35 0.06 0.35      "unlist" 0.06 0.35 0.06 0.35      "print.default" 0.02 0.12 0.02 0.12      "print" 0.02 0.12 0.00 0.00            $sample.interval      [1] 0.02            $sampling.time      [1] 17.1                  C'est donc cbind() qui consomme le plus de temps machine et donc powers1() par voie de conséquence... Chacun a en général ses habitudes de programmation, de nommage et d'indentation, mais il semble qu'il y ait des consensus, sachant que le principal critère d'écriture de code doit être la lisibilité. S'il est évident qu'il ne faut pas utiliser des noms de variables à une seule lettre, des noms trop longs ne facilitent en rien la lecture. On a principalement le choix entre le style chameau, le style tout collé, le style avec des points et le style tout_souligné comme par exemple : # différentes façons de nommer la variable      # "nombre de lignes" du dataframe mesData            n <- nrow(mesData) # à la fainéante (1)      nb <- nrow(mesData) # à la fainéante (2)      nl <- nrow(mesData) # à la fainéante (3)            nbl <- nrow(mesData) # en mode court (3 lettres)      nlig <- nrow(mesData) # en mode court (4 lettres)            nblig <- nrow(mesData) # en mode tout collé      nbLig <- nrow(mesData) # en mode chameau (camelCase)      nb.lig <- nrow(mesData) # en mode avec des points      nb_lig <- nrow(mesData) # en mode souligné            nombre.de.lignes <- nrow(mesData) # en mode avec des points très long      nombredelignes <- nrow(mesData) # en mode tout collé très long      nombreDeLignes <- nrow(mesData) # en mode chameau (camelCase) très long      nombre_de_lignes <- nrow(mesData) # en mode souligné très long                  Dans la mesure où R utilise des points pour la notation objets, il vaut peut-être mieux éviter le style avec des points pour nommer les fonctions, mais pas pour les variables... La notation hongroise a eu ses adeptes pendant un moment et peut permettre de distinguer rapidement les variables numériques et leur équivalent texte formaté, par exemple. # exemple de notation hongroise adaptée :      # n_ pour entier et f_ pour sortie formatée            n_lig <- nrow( mesData)            f_lig <- sprintf("%03d",n_lig)            Il y a des pages Web qui référencent ou conseillent des pratiques de style pour R. On peut citer notamment la recommandation Google, la RCC, celle de stat405. Une pratique classique peut être d'utiliser des variables à trois lettres à condition de founir dès le début de la fonction une table de correspondance : # exemple de table de correspondance            analyseMalah <- function( dfMalah,...) {            ###############################################            # nbl : nombre de lignes      # nbc : nombre de colonnes      # nbv : nombre de variables      # nbm : nombre de malades            # vpr : variable principale      # idpr : son indice      # rpr : son R2      ...                  Enfin, il ne faut pas confondre saisie et codage. Rien n'interdit d'utiliser des variables i, j, n... sous éditeur, de mettre au point la fonction puis de faire des édition/remplacer pour obtenir des noms de variables explicites : ##############################################################      #      # en cours de développement      #      #############################################################            concatv <- function(ldtc) {                  lg <- length(ldtc)      mm <- max(sapply(X=ldtc,FUN=length))      mr <- matrix("",nrow=mm,ncol=lg)            for (i in seq(ldtc)) {      el <- ldtc[[i]]      le <- length(el)      mr[ (1:le), i ] <- el      } # fin pour i            print(mres,quote=FALSE)            } # fin de fonction concatv            ##############################################################      #      # une fois la fonction mise au point      #      ##############################################################            concatv <- function(ldtc) {                  lngList <- length(ldtc)      lngMax <- max(sapply(X=ldtc,FUN=length))      mres <- matrix("",nrow=lngMax,ncol=lngList)            for (indElt in seq(ldtc)) {      elt <- ldtc[[indElt]]      lngElt <- length(elt)      mres[ (1:lngElt), indElt ] <- elt      } # fin pour indElt            print(mres,quote=FALSE)            } # fin de fonction concatv            Rappelons quand même quelques précautions élémentaires de codage : écrire les commentaires avant les instructions ; mettre des accolades même pour un corps à une seule instruction ; commenter ce que fait chaque fonction ; prévoir les erreurs usuelles ; penser aux structures vides ; mettre des parenthèses pour lever toute ambiguité... Il est en général conseillé d'indenter avec au moins trois espaces les structures, c'est-à-dire de laisser au moins trois espaces à droite pour les boucles, les tests, les fonctions... Nous recommandons aussi de nommer les fins de structures au lieu de laisser les accolades anonymes : même correcte, une écriture comme }}} }} est toujours peu lisible. Enfin, une écriture comme # fin de fonction maFonc permet de repérer informatiquement le corps complet de la fonction, plutôt que d'essayer de retrouver la parenthèse fermante correspondant à maFonc <- function(....) {. Nous conseillons d'écrire en fin fonction plutôt que function car cela permet rapidement de repérer en ligne de commande avec grep les débuts de fonctions... Cela permet aussi facilement de mettre en oeuvre une page web avec code-source de la fonction et des exemples d'utilisation comme pour statgh.r. # mauvais style : pas d'indendation            findruns <- function(x,k) {      n <- length(x)      runs <- NULL      for (i in 1:(n-k+1)) {      if (all(x[i:(i+k-1)]==1)) {      runs <- c(runs,i)      }} return(runs) }            # un style meilleur : indendation      # et nommage des fins de structure,      # texte plus aéré...            ##################################################            findruns <- function(x,k) {            ##################################################            n <- length(x)      runs <- NULL            for (i in 1:(n-k+1)) {      if (all(x[i:(i+k-1)]==1)) {      runs <- c(runs,i)      } # fin si      } # fin pour i            return(runs)            } # fin de fonction findruns            ##################################################                  C'est si simple de détecter les fonctions et leurs exemples avec ce simple marquage que nous fournissons une "boite vide" à l'adresse  rfns  pour y déposer vos fonctions et vos exemples. Lorsqu'on développe principalement pour soi, on peut se fixer des standards d'écriture, comme par exemple i désignera toujours un indice de vecteur, ilig un indice de ligne, jcol un indice de colonne, pdv une plage de valeurs... Il est toujours bon de fournir de l'aide, le rappel de la syntaxe, les cas particuliers, les objets retournés, les précautions d'emploi... Fournir des jeux de données d'exemples (surtout si on utilise des structures de données un peu techniques) est bien sûr recommandé, de même que des exemples d'utilisations et de sorties... Ce genre d'aide est utile même quand on développe pour soi : il arrive que l'on reprenne un programme six mois plus (quand ce n'est pas un an...) et retrouver ce qu'on avait fait, quelle logique on avait suivie, quels choix avaient été programmés sont plus faciles à relire quand ils ont été rédigés qu'à retrouver en parcourant toutes les lignes de code... Quand on écrit du code, que ce soit en R ou dans un autre langage, il arrive qu'on ait à modifier le code, à le réutiliser, le réarranger. Avec R, on vient souvent ajouter des nouveaux paramètres. Si la fonction doit réaliser des tracés graphiques, une bonne habitude à prendre est d'ajouter l'ellipse notée ... ; réviser l'exercice 8 de la session 1 pour plus de détails. La PO (programmation objet) et la POO (programmation orientée objet) sont des styles de programmation adaptés au développement d'importance. Avec des objets, la programmation est en principe "plus propre" et plus facile à maintenir. Le concept de base est de mettre ensemble, de façon protégée et cohérente les données et les fonctions liées aux mêmes structures de données. Les objets de type S3 sont issus du langage S et sont les plus anciens. Les "nouveaux" objets sont de type S4. Beaucoup d'objets statistiques implémentés dès le départ (comme les résultats de régression linéaire) sont de type S3. La définition et la gestion des objets de type S3 et S4 sont différentes et relativement compatibles... On pourra lire avec profit le petit manuel d'introduction aux objets S4 en R nommé A (Not So) Short Introduction to S4 et dont la traduction française est Petit Manuel de S4, documents disponibles dans la page des contributions de documentation sur le site du CRAN à la rubrique des manuels. On y trouve le vocabulaire classique de classe, attribut, méthode, héritage, encapsulation, polymorphisme, La première chose à faire est de bien réfléchir à l'ensemble des fonctionnalités que l'on veut implémenter, quitte à tout prévoir mais à ne pas tout implémenter tout de suite. Comme des classes peuvent contenir des sous-classes, il faut avec prudence et clarté définir ce qui appartient à chaque classe, en terme de données (attributs) ou en terme de fonctions (méthodes) sachant que le polymorphisme permet d'utiliser le même nom pour des fonctions différentes (car dans différentes classes) réalisant des actions similaires. Pour notre problème, la classe des VS (variables statistiques) contiendra la sous-classe des VQT (variables statistiques quantitatives) et la sous-classe des VQL (variables statistiques qualitatives). Nous allons supposer qu'en toute généralité, une VS est définie par un titre, une référence et des valeurs numériques. Pour une VQT, il faut adjoindre l'unité et pour une VQL il faut ajouter les modalités. En voici une pseudo-définition textuelle où on ne détaille pas pour l'instant les paramètres des fonctions : +--- Classe VS (variable statistique)            |      |      | Attribut titre chaine_de_caractères # titre général      | Attribut ref chaine_de_caractères # référence      | Attribut valeurs vecteur_numérique # "les données"      |      | Méthode print # pour l'affichage      | Méthode dsc # dsc pour "décrit"      |      |            -------+ -- Sous-Classe VQT (variable statistique quantitative)            |      |      | Attribut unite chaine_de_caractères      |      | Méthode print # pour l'affichage      | Méthode plot # pour le tracé      | Méthode dsc # dsc pour "décrit"      |      |            -------+ -- Sous-Classe VQT (variable statistique qualitative)            |      |      | Attribut modalites vecteur_de_chaines_de_caractères      | Attribut facteur # au sens de R      |      | Méthode print # pour l'affichage      | Méthode plot # pour le tracé      | Méthode dsc # dsc pour "décrit"      | Méthode as.vector # pour l'export      |      |      ...            On définit les classes en R via setClass() et on instancie les objets via new() sachant que toutes les fonctions pour les objets sont regroupées dans le package nommé methods. Avec la définition d'objets, R permet de surveiller les types de données utilisés, via la notion de representation(). Dans la mesure où les classes S4 sont plus "évoluées" que les classes S3 et sont aujourd'hui (2013) recommandées, nous donnons la définition en S4 : ## objets S4, classe VS (variables statistiques)            # 1. définition            setClass(            Class="VS",            representation=representation(      titre="character",      ref="character",      valeurs="numeric"      )            ) # fin de définition de la classe VS            # 2. utilisation            # refusé :      # badVar1a <- new("VS",titre="ages",valeurs="30 40 50 60")      # non refusé et pourtant il manque ref :      # badVar1b <- new("VS",titre="ages",valeurs=c(30,40,50,60))            # accepté :            var1 <- new(Class="VS",titre="ages",ref="gh DEMO",valeurs=c(30,40,50,NA,60,30))            # ce qu'on peut en connaitre            print( class(var1) )      print( typeof(var1) )      print( attributes(var1) )      print( slotNames(var1) )      print( getSlots( class(var1)) )      print( getClass( class(var1)) )            print( var1@titre )      print( attr(var1,"titre") )      # not run : print(var1$titre)                        > ## objets S4      >      > # 1. définition      >      > setClass(      +      + Class="VS",      +      + representation=representation(      + titre="character",      + ref="charac ..." ... [TRUNCATED]      [1] "VS"            > # 2. utilisation      >      > # refusé :      > # badVar1a <- new("VS",titre="ages",valeurs="30 40 50 60")      > # Erreur dans validObject(.Object) :      > # objet de classe "VS" incorrect :      > # invalid object for slot "valeurs" in class "VS":      > # got class "character", should be or extend class "numeric"      > #      > # non refusé et pourtant il manque ref :      > # badVar1b <- new("VS",titre="ages",valeurs=c(30,40,50,60))            > # ce qu'on peut en connaitre      >      > print( class(var1) )      [1] "VS"      attr(,"package")      [1] ".GlobalEnv"            > print( typeof(var1) )      [1] "S4"            > print( attributes(var1) )      $titre      [1] "ages"            $ref      [1] "gh DEMO"            $valeurs      [1] 30 40 50 NA 60 30            $class      [1] "VS"      attr(,"package")      [1] ".GlobalEnv"                  > print( slotNames(var1) )      [1] "titre" "ref" "valeurs"            > print( getSlots( class(var1)) )      titre ref valeurs      "character" "character" "numeric"            > print( getClass( class(var1)) )      Class "VS" [in ".GlobalEnv"]            Slots:            Name: titre ref valeurs      Class: character character numeric            > print( var1@titre )      [1] "ages"            > print( attr(var1,"titre ") )      [1] "ages"            > # not run : print(var1$titre)      > # Erreur dans var1$titre : $ operator not defined for this S4 class            Le vocabulaire de R utilise la notion de slot pour gérer les attributs. La notation @ permet d'accéder aux informations associées. R, via la notion de representation(), permet de définir la signature et le typage des données de l'objet. A l'aide de prototype et validity, on peut garantir à la création que les données de l'objet sont bien typées, et qu'elles sont fonctionnelles (non vides, plus grandes que, etc.). Voici donc une implémentation plus complète : ## objets S4, définition plus complète pour la classe VS (variables statistiques)            setClass(            Class="VSv2",            representation=representation(      titre="character",      ref="character",      valeurs="numeric"      )      ,      prototype=prototype(      titre="",      ref="",      valeurs=numeric(0)      )      ,      validity=function(object) {      (nchar(object@titre)>0) &&      (nchar(object@ref)>0)      } # fin de validity            ) # fin de définition de la classe VSv2            ## utilisation            # refusé désormais :      # badVar2 <- new("VSv2",titre="ages",valeurs=c(30,40,50,60))            # accepté :            var2 <- new(Class="VSv2",titre="ages",ref="gh DEMO",valeurs=c(30,40,50,NA,60,30))                  > ## objets S4, définition plus complète      >      > setClass(      +      + Class="VSv2",      +      + representation=representation(      + titre="character",      + ref= .... [TRUNCATED]      [1] "VSv2"            > ## utilisation      >      > # refusé désormais :      > # badVar2 <- new("VSv2",titre="ages",valeurs=c(30,40,50,60))      > # Erreur dans validObject(.Object) : objet de classe "VSv2" incorrect : FALSE            > var2 <- new(Class="VSv2",titre="ages",ref="gh DEMO",valeurs=c(30,40,50,NA,60,30))                  Mais il est possible de définition de façon extérieure à la création la vérification de la validité d'un objet avec setValidity() et validObject(). De plus, pour afficher un objet, on peut redéfinir la fonction show() pour la classe considérée  : ## objets S4 avec validation plus explicite et gestion de l'affichage            setClass(      Class="VSv3",      representation=representation(      titre="character",      ref="character",      valeurs="numeric"      ) ,      prototype=prototype(      titre="",      ref="",      valeurs=numeric(0)      ) ,      validity=function(object) {      (nchar(object@titre)>0) &&      (nchar(object@ref)>0)      } # fin de validity            ) # fin de définition de la classe VSv3            setValidity(      Class="VSv3",      method= function(object) {      ok <- TRUE      if (nchar(object@titre)==0) {      cat("Erreur : il manque le titre \n")      ok <- FALSE      } # finsi      if (nchar(object@ref)==0) {      cat("Erreur : il manque la référence \n")      ok <- FALSE      } # finsi      return( ok )      } # fin de méthode      ) # fin de fonction de validation pour VSv3            setMethod(      signature="VSv3",      f="show",      definition=function(object) {      cat("Variable statistique : ",object@titre," références : ",object@ref,"\n")      nbval <- length(object@valeurs)      nbvaldis <- length(unique(object@valeurs))      cat(" avec",nbval,"valeurs")      if (nbval==nbvaldis) {      cat(" uniques")      } else {      cat(" dont ",nbvaldis,"uniques")      } # fin si      nbna <- sum(is.na(object@valeurs))      if (nbna>0) {      cat(" et",nbna,"valeur(s) manquante(s)")      } # fin si      cat(".\n")      } # fin de définition      ) # fin de setMethod show pour VSv3            var3a <- new(Class="VSv3",titre="ages",ref="gh DEMO",valeurs=c(30,40,50,60))      show(var3a)      print(var3a)            var3b <- new(Class="VSv3",titre="ages",ref="gh DEMO",valeurs=c(30,40,50,NA,60,30,60,30))      show(var3b)      print(var3b)            # not run :      # badVar3a <- new("VSv3",titre="ages",valeurs="30 40 50 60")      # badVar3b <- new("VSv3",titre="ages",valeurs=c(30,40,50,60))                  > ## objets S4 avec validation plus explicite et gestion de l'affichage      >      > setClass(      + Class="VSv3",      + representation=representation(      + titr .... [TRUNCATED]      [1] "VSv3"            > setValidity(      + Class="VSv3",      + method= function(object) {      + ok <- TRUE      + if (nchar(object@titre)==0) {      + cat("Erreur : il manqu ..." ... [TRUNCATED]      Class "VSv3" [in ".GlobalEnv"]            Slots:            Name: titre ref valeurs      Class: character character numeric            > setMethod(      + signature="VSv3",      + f="show",      + definition=function(object) {      + cat("Variable statistique : ",object@titre," références : ",o .... [TRUNCATED]      [1] "show"            > var3a <- new(Class="VSv3",titre="ages",ref="gh DEMO",valeurs=c(30,40,50,60))            > show(var3a)      Variable statistique : ages références : gh DEMO      avec 4 valeurs uniques.            > print(var3a)      Variable statistique : ages références : gh DEMO      avec 4 valeurs uniques.            > var3b <- new(Class="VSv3",titre="ages",ref="gh DEMO",valeurs=c(30,40,50,NA,60,30,60,30))            > show(var3b)      Variable statistique : ages références : gh DEMO      avec 8 valeurs dont 5 uniques et 1 valeur(s) manquante(s).            > print(var3b)      Variable statistique : ages références : gh DEMO      avec 8 valeurs dont 5 uniques et 1 valeur(s) manquante(s).            > badVar3a <- new("VSv3",titre="ages",valeurs="30 40 50 60")      Erreur dans validObject(.Object) :      invalid class "VSv3" object: invalid object for slot "valeurs" in class "VSv3":      got class "character", should be or extend class "numeric"            > badVar3b <- new("VSv3",titre="ages",valeurs=c(30,40,50,60))      Erreur : il manque la référence      Erreur dans validObject(.Object) : invalid class "VSv3" object: FALSE            Pour définir des sous-classes, il suffit d'utiliser contains dans la définition de la classe. Les sous-classes héritent les attributs et les méthodes de la sur-classe, il ne faut donc définir que ce qui est ajouté. ## objets S4 : sous-classes VSqt (quantitatives) et VSql (qualitatives)      ## de la classe VS (variables statistiques)            setClass(            Class="VSqt",      contains="VSv3",            representation=representation(      unite="character"      )      ,      prototype=prototype(      unite="???"      )            ) # fin de définition de la classe VSqt            setClass(            Class="VSql",            representation=representation(      modalites="character",      donnees="factor"      )      ,      contains="VSv3"            ) # fin de définition de la classe VSqt            var4qt <- new(Class="VSqt",titre="ages",ref="gh DEMO",valeurs=c(30,40,50,NA,60,30),unite="ans")            print(var4qt)            print( slotNames(var4qt) )      print( getSlots( class(var4qt)) )      print( getClass( class(var4qt)) )                        > ## objets S4 : sous-classes      >      > setClass(      +      + Class="VSqt",      + contains="VSv3",      +      + representation=representation(      + unite="character"      + .... [TRUNCATED]      [1] "VSqt"            > setClass(      +      + Class="ghql",      +      + representation=representation(      + modalites="character"      + )      + ,      + contains="VSv3"      +      + ) # fin de déf .... [TRUNCATED]      [1] "ghql"            > var4qt <- new(Class="VSqt",titre="ages",ref="gh DEMO",valeurs=c(30,40,50,NA,60,30),unite="ans")            > print(var4qt)      Variable statistique : ages références : gh DEMO      avec 6 valeurs dont 5 uniques et 1 valeur(s) manquante(s).            > print( slotNames(var4qt) )      [1] "unite" "titre" "ref" "valeurs"            > print( getSlots( class(var4qt)) )      unite titre ref valeurs      "character" "character" "character" "numeric"            > print( getClass( class(var4qt)) )      Class "VSqt" [in ".GlobalEnv"]            Slots:            Name: unite titre ref valeurs      Class: character character character numeric            Extends: "VSv3"            Avec la fonction initialize on peut effectuer tous les traitements nécessaires à la bonne initialisation de l'objet.       # initialisateur pour la classe ghql, on crée le facteur                  setMethod(      signature="ghql",      f="initialize",      definition=function(object) {            ...            # forcer à tester la validation      # ne pas oublier de renvoyer l'objet            } # fin de définition            ) # fin de setMethod                  Enfin, il est usuel de définir une fonction qui porte le nom de la classe et qui initialise «ce qu'il faut où il faut» : # ajout d'une fonction du même nom que la classe            VSqt <- function(sonTitre,saRef,sesValeurs=numeric(0),sonUnite="???") {            if (missing(sonTitre)) {      stop("pour créer un objet VSqt, il faut impérativement donner son titre\n")      } # finsi            if (missing(saRef)) {      stop("pour créer un objet VSqt, il faut impérativement donner sa référence\n")      return()      } # finsi            return( new(Class="VSqt",titre=sonTitre,ref=saRef,valeurs=sesValeurs,unite=sonUnite) )            } # finsi            # à tester :            v6a <- VSqt(sonTitre="poids",saRef="dossier HER (gh)",sesValeurs=c(8,12,15),sonUnite="kg")      print( v6a )            v6b <- VSqt("taille"," exemple GH")      print( v6b )            v6c <- VSqt()                  > # ajout d'une fonction du même nom que la classe      >      > VSqt <- function(sonTitre,saRef,sesValeurs=numeric(0),sonUnite="???") {      +      + if (missing(s .... [TRUNCATED]            > # à tester :      >      > v6a <- VSqt(sonTitre="poids",saRef="dossier HER (gh)",sesValeurs=c(8,12,15),sonUnite="kg")            > print( v6a )      Variable statistique : poids références : dossier HER (gh)      avec 3 valeurs uniques.            > v6b <- VSqt("taille"," exemple GH")            > print( v6b )      Variable statistique : taille références : exemple GH      avec 0 valeurs uniques.            > v6c <- VSqt()      Erreur dans VSqt() :      pour créer un objet VSqt, il faut impérativement donner son titre                  Pour accéder aux valeurs, que ce soit en lecture ou en écriture, on écrit en général des fonctions get et set et R permet aussi de définir des fonctions <- : ##########################################      #      # fonctions pour accéder aux slots      #      ##########################################            setGeneric("getTitre",      function(object) { standardGeneric("getTitre") }      ) # fin de setGeneric pour getTitre            setGeneric("getRef",      function(object) { standardGeneric("getRef") }      ) # fin de setGeneric pour getTitre            setGeneric("getValeurs",      function(object) { standardGeneric("getValeurs") }      ) # fin de setGeneric pour getValeurs            setGeneric("getUnite",      function(object) { standardGeneric("getUnite") }      ) # fin de setGeneric pour getUnite            # ------------------------------------------------            setMethod("getTitre","VSqt",      function(object) { return( object@titre ) }      ) # fin de setMethod pour getTitre            setMethod("getRef","VSqt",      function(object) { return( object@ref ) }      ) # fin de setMethod pour getRef            setMethod("getValeurs","VSqt",      function(object) { return( object@valeurs ) }      ) # fin de setMethod pour getValeurs            setMethod("getUnite","VSqt",      function(object) { return( object@unite ) }      ) # fin de setMethod pour getUnite            # ------------------------------------------------            setGeneric("setTitre<-",      function(object,...) { standardGeneric("setTitre<-") }      ) # fin de setGeneric pour setTitre            setGeneric("setRef<-",      function(object,...) { standardGeneric("setRef<-") }      ) # fin de setGeneric pour setTitre            setGeneric("setValeurs<-",      function(object,...) { standardGeneric("setValeurs<-") }      ) # fin de setGeneric pour setValeurs            setGeneric("setUnite<-",      function(object,...) { standardGeneric("setUnite<-") }      ) # fin de setGeneric pour setUnite            # ------------------------------------------------            setReplaceMethod(      f="setTitre",      signature="VSqt",      definition=function(object,value) {      object@titre <- value      return(object)      }      ) # fin de setReplaceMethod pour setTitre            setReplaceMethod(      f="setRef",      signature="VSqt",      definition=function(object,value) {      object@ref <- value      return(object)      }      ) # fin de setReplaceMethod pour setRef            setReplaceMethod(      f="setValeurs",      signature="VSqt",      definition=function(object,value) {      object@valeurs <- value      return(object)      }      ) # fin de setReplaceMethod pour setValeurs            setReplaceMethod(      f="setUnite",      signature="VSqt",      definition=function(object,value) {      object@unite <- value      return(object)      }      ) # fin de setReplaceMethod pour setUnite            ##########################################      #      # exemples :      #      ##########################################            v7 <- VSqt(sonTitre="poids",saRef="dossier HER (gh)",sesValeurs=c(8,12,15),sonUnite="kg")      print( v7 )            getTitre( v7 )      getRef( v7 )      getValeurs( v7 )      getUnite( v7 )            setTitre( v7 ) <- "poidsG"      setRef( v7 ) <- "transf. en grammes de pois HER "      setValeurs( v7 ) <- getValeurs(v7)*1000      setUnite( v7 ) <- "g"            print( v7 )      getValeurs( v7 )                        > ##########################################      > #      > # fonctions pour accéder aux slots      > #      > ##########################################      >      .... [TRUNCATED]            > ##########################################      > #      > # exemples :      > #      > ##########################################      >      > v7 <- VSqt(sonTitre="poids",saRe .... [TRUNCATED]            > print( v7 )      Variable statistique : poids références : dossier HER (gh)      avec 3 valeurs uniques.            > getTitre( v7 )      [1] "poids"            > getRef( v7 )      [1] "dossier HER (gh)"            > getValeurs( v7 )      [1] 8 12 15            > getUnite( v7 )      [1] "kg"            > setTitre( v7 ) <- "poidsG"      > setRef( v7 ) <- "transf. en grammes de pois HER "      > setValeurs( v7 ) <- getValeurs(v7)*1000      > setUnite( v7 ) <- "g"            > print( v7 )      Variable statistique : poidsG références : transf. en grammes de pois HER      avec 3 valeurs uniques.            > getValeurs( v7 )      [1] 8000 12000 15000            Au-delà de tout ce travail, le grand avantage de la programmation objet est de fournir une grande souplesse d'utilisation parce que tout est masqué et cohérent (quand c'est bien programmé) : # traitement objet d'une qt            print( maqt )      plot( maqt )      decrit( maqt )            # traitement objet d'une ql            print( maql )      plot( maql )      decrit( maql )            Pour finir ce petit tour des objets S4, une implémentation complète de classes d'objets se doit aussi de fournir des fonctions de test de classe comme is.VS et des fonctions de conversion de classe comme as.VS... Ce n'est pas très difficile de produire un package, une fois qu'on a écrit ses fonctions, préparé des jeux de données d'essai et qu'on a pensé aux exemples à fournir. R dispose d'une fonction package.skeleton() qui prépare tout le travail pour nous. Nous présentons ici une version simple de package à réaliser qui n'utilise pas de programme en C ou dans un autre langage qui viendrait s'ajouter aux fonctions et aux données R. C'est donc un exemple minimaliste et simpliste. Admettons par exemple qu'on veuille inventer la fonction cdv (coefficient de variation) et utiliser comme jeu d'essai les données du fichier elfdemo.dat pour créer le package demoGh sous Windows, sur le disque Z:  : cdv <- function(x) sd(x)/mean(x)      elf <- read.table("elfdemo.dat",head=TRUE,row.names=1)      package.skeleton(name="demoGh",path="Z:/",list=c("cdv","elf"))            Création des répertoires...      Création de DESCRIPTION...      Création de 'Read-and-delete-me'...      Sauvegarde des fonctions et des données...      Création des fichiers d'aide ...      Terminé.      Les étapes suivantes sont décrites dans 'Z://demoGh/Read-and-delete-me'.            R vient alors créer le dossier Z://demoGh/ et y crée toutes les entrées à compléter. On trouvera dans l'archive demoGh.zip les dossiers et les fichiers créés... On pourra en particulier consulter demoGh-package.Rd.txt, elf.Rd.txt, cdv.Rd.txt et suivre les instructions de Read-and-delete-me.txt . Avec beaucoup d'impatience, on peut essayer tout de suite de construire le package en ligne de commande avec Z:\>"C:\Program Files\R\R-2.13.2\bin\i386\Rcmd.exe" check demoGh            Rcmd nous prévient alors que certains champs sont vides et d'autres mal remplis : * installing *source* package 'demoGh' ...      ** R      ** data      ** preparing package for lazy loading      ** help      Avis : Z:/demoGh/man/demoGh-package.Rd:34: All text must be in a section      Avis : Z:/demoGh/man/demoGh-package.Rd:35: All text must be in a section      *** installing help indices      Erreur dans Rd_info(db[[i]]) :      champ \title manquant ou vide dans 'Z:/demoGh/man/cdv.Rd'      les fichiers Rd doivent avoir un \title non vide                  Et il faut alors aller consulter le fichier Z:\demoGh.Rcheck\00install.out ou son équivalent pour Linux... ainsi que le détail de la vérification dans le fichier 00check.log qui correspond à ce qui suit : * using log directory 'Z://demoGh.Rcheck' * using R version 2.13.2 (2011-09-30) * using platform: i386-pc-mingw32 (32-bit) * using session charset: ISO8859-1 * checking for file 'demoGh/DESCRIPTION' ... OK * checking extension type ... Package * this is package 'demoGh' version '1.0' * checking package dependencies ... OK * checking if this is a source package ... OK * checking for .dll and .exe files ... OK * checking whether package 'demoGh' can be installed ... OK * checking package directory ... OK * checking for portable file names ... OK * checking DESCRIPTION meta-information ... WARNING Non-standard license specification: What license is it under? Standardizable: FALSE * checking top-level files ... OK * checking index information ... OK * checking package subdirectories ... OK * checking R files for non-ASCII characters ... OK * checking R files for syntax errors ... OK * checking whether the package can be loaded ... OK * checking whether the package can be loaded with stated dependencies ... OK * checking whether the package can be unloaded cleanly ... OK * checking for unstated dependencies in R code ... OK * checking S3 generic/method consistency ... OK * checking replacement functions ... OK * checking foreign function calls ... OK * checking R code for possible problems ... OK * checking Rd files ... WARNING Error : demogh-package.rd: non-ASCII input and no declared encoding problem found in 'demogh-package.rd' * checking Rd metadata ... OK * checking Rd cross-references ... OK * checking for missing documentation entries ... OK * checking for code/documentation mismatches ... OK * checking Rd \usage sections ... OK * checking Rd contents ... OK * checking for unstated dependencies in examples ... OK * checking contents of 'data' directory ... OK * checking data for non-ASCII characters ... OK * checking data for ASCII and uncompressed saves ... OK * checking examples ... OK * checking PDF version of manual ... OK WARNING: There were 2 warnings, see 'Z://demoGh.Rcheck/00check.log' for details Après avoir lu http://cran.r-project.org/doc/manuals/R-exts.html ou son équivalent PDF, on comprend enfin(!) qu'il faut remplir les champs marqués ~ et % des fichiers .rd (Windows) ou .Rd (Linux). Voici les trois fichiers, rapidement complétés : Fichier man/demoGh-package.rd modifié, l'original est ici : \name{demoGh-package} \alias{demoGh-package} \alias{demoGh} \docType{package} \title{ exemple pour l'Institut Pasteur (décembre 2013) } \description{ contient juste la fonction cdv et le jeu de données ELF. } \details{ \tabular{ll}{ Package: \tab demoGh\cr Type: \tab Package\cr Version: \tab 1.0\cr Date: \tab 2012-03-08\cr License: GNU LazyLoad: \tab yes\cr } ~~ An overview of how to use the package, including the most important ~~ ~~ functions ~~ } \author{ (gH) Maintainer: Who to complain to <yourfault@somewhere.net> (gH) } \references{ (gH) } ~~ Optionally other standard keywords, one per line, from file KEYWORDS in ~~ ~~ the R documentation directory ~~ \keyword{ package } \seealso{ mean } \examples{ data(elf) cdv ( elf$AGE ) data(cars) cdv(cars$dist) } Fichier man/elf.rd modifié, l'original est ici : \name{elf} \alias{elf} \docType{data} \title{ELF datafile } \description{ French dataset (1991) of 99 people with gender and education status. } \usage{data(elf)} \format{ A data frame with 99 observations on the following 3 variables. \describe{ \item{\code{SEXE}}{a numeric vector} \item{\code{AGE}}{a numeric vector} \item{\code{ETUD}}{a numeric vector} } } \details{ 0 means male and 1 female for the SEXE variable. AGE refers to the age of the person, in years. ETUD is a code for education level... } \source{ (gH) } \references{ (gH) } \examples{ data(elf) cdv( elf$AGE ) } \keyword{datasets} Fichier man/cdv.rd modifié, l'original est ici : \name{cdv} \alias{cdv} %- Also NEED an '\alias' for EACH other topic documented here. \title{ Computes the coefficient of variation, that is, the ratio sd(x)/mean(x). %% ~~function to do ... ~~ } \description{ coefficient de variation } \usage{ cdv(x) } %- maybe also 'usage' for other objects documented here. \arguments{ \item{x}{ %% ~~Describe \code{x} here~~ } } \details{ calcul minimal } \value{ (gH) } \references{ (gH) } \author{ (gH) } \note{ 18/20 ! } %% ~Make other sections like Warning with \section{Warning }{....} ~ \seealso{ \code{\link{help}} } \examples{ data(elf) cdv( elf$AGE ) data(cars) cdv(cars$dist) } % Add one or more standard keywords, see file 'KEYWORDS' in the % R documentation directory. \keyword{ ~mean } \keyword{ ~sd } Le package est alors presque terminé (hum !) et R arrive à produire le fichier PDF d'aide, les fichiers d'aide... Z:\>"C:\Program Files\R\R-2.13.2\bin\i386\Rcmd.exe" check demoGh      Z:\>"C:\Program Files\R\R-2.13.2\bin\i386\Rcmd.exe" build demoGh            Consultez l'archive demoGh2.zip pour voir tout ce qui a été rajouté, comme demoGh-manual.pdf et le fichier 00install.out.txt dont le contenu suit * installing *source* package 'demoGh' ... ** R ** data ** preparing package for lazy loading ** help Avis : Z:/demoGh/man/demogh-package.rd:32: All text must be in a section Avis : Z:/demoGh/man/demogh-package.rd:33: All text must be in a section *** installing help indices ** building package indices ... ** testing if installed package can be loaded * DONE (demoGh) Les commandes R CMD check demoGh et R CMD build demoGh, lorsque le check ne renvoie plus d'erreurs ni de warnings, montrent qu'il est relativement simple du point de vue technique de construire un package. C'est par contre un travail long, soigné, pour lequel il faut être trés vigilant, rigoureux et précis. Cela oblige à documenter chaque fonction, chaque jeu de données, à fournir des exemples, des mots-clés, des renvois, des références... Il existe d'autres façons de faire, comme de mettre du code en ligne, des codes-sources et des exemples à télécharger... C'est ce que nous proposons avec les fonctions du fichier statgh.r complété par un fichier d'exemples d'utilisation et une interface Web pour consulter les fonctions. Un package est la forme la plus aboutie mais aussi la plus figée des fonctions et des données. Pour approfondir le sujet, une documentation en français est celle de C. Genolini. Chaque famille de fonction (tests, régression...) a sa propre classe d'objets, comme on peut le voir ci-dessous.       # données            library(MASS)            data(cars)      data(iris)      data(iris3)            pg <- lit.dar("pg.dar")      iris4 <- iris[,-5]      Iris <- data.frame(rbind(iris3[,,1], iris3[,,2], iris3[,,3]),      Sp = rep(c("s","c","v"), rep(50,3)))      train <- sample(1:150, 75)            # objets issus de tests            tt <- t.test(cars$dist)      tw <- wilcox.test(cars$dist)      tc <- cor.test(cars$dist,cars$speed)            # objets issus de régression            rlis <- with(cars,lm(speed~dist))      rlim <- lm(Petal.Width~.,data=iris4)      rlob <- with(pg,glm(GROUPE~TAILLE,family="binomial"))      rlda <- lda(Sp ~ ., Iris, prior = c(1,1,1)/3, subset = train)            # autres objets            hdd <- hist(cars$dist,plot=FALSE)      arl <- anova(rlis)      sdf <- summary(iris)      tap <- table(iris$Species)            # analyse de la liste des objets            ldo <- list(tt,tw,tc,rlis,rlim,rlob,rlda,hdd,arl,sdf,tap)      nbo <- length(ldo)            mres <- matrix(nrow=nbo,ncol=3)      mres <- data.frame(mres)            row.names(mres) <- c(      "Test t",      "Test de Wilcox",      "Test de corrélation",      "Régression linéaire simple",      "Régression linéaire multiple",      "Régression logistique binaire",      "Analyse discriminante linéaire",      "Histogramme",      "Anova de régression linéaire",      "Summary de data frame",      "Tri à plat"      ) # fin de row.names            colnames(mres) <- c(      "Classe(s)",      "Liste?",      "Longueur"      ) # fin de colnames            mres[,1] <- unlist(sapply(X=ldo,FUN=function(x) paste(class(x),collapse=" ")))      mres[,2] <- unlist(sapply(X=ldo,FUN=is.list))      mres[,3] <- unlist(sapply(X=ldo,FUN=length))            mres[,1] <- reformate(mres[,1])            print(mres,quote=FALSE)            Classe(s) Liste? Longueur      Test t htest __________ TRUE 9      Test de Wilcox htest __________ TRUE 7      Test de corrélation htest __________ TRUE 9      Régression linéaire simple lm _____________ TRUE 12      Régression linéaire multiple lm _____________ TRUE 12      Régression logistique binaire glm lm _________ TRUE 30      Analyse discriminante linéaire lda ____________ TRUE 10      Histogramme histogram ______ TRUE 6      Anova de régression linéaire anova data.frame TRUE 5      Summary de data frame table __________ FALSE 30      Tri à plat table __________ FALSE 3            Ce n'est certainement pas une bonne idée que de vouloir comprendre des calculs statistiques en lisant du code informatique qui les implémente. Mais en admettant qu'on veuille juste reproduire le calcul des probabilités a posteriori produit par predict(lda(...),type="class"), il suffit de récupérer le code-source de la fonction sous-jacente. Comme la fonction lda() du package MASS renvoie des objets de type lda, c'est donc la fonction predict.lda() qu'il faut utiliser. Malheureusement, elle est "invisible" (ou plus exactement masquée) et on ne peut donc la voir ni dans les aides ni dans la session courante de R : > apropos("predict")            "makepredictcall" "napredict" "predict" "predict.glm" "predict.lm"            > library(MASS)            > ls(name="package:MASS",pattern="lda")            "lda" "ldahist"            > ls(name="package:MASS",pattern="predict")            character(0)            Pour obtenir le code-source de la fonction predict.lda(), il faut utiliser getAnywhere() du package utils : > help(getAnywhere)            getAnywhere package:utils R Documentation            Retrieve an R Object, Including from a Namespace            Description:            These functions locate all objects with name matching their      argument, whether visible on the search path, registered as an S3      method or in a namespace but not exported. \u2018getAnywhere()'      returns the objects and \u2018argsAnywhere()' returns the arguments of      any objects that are functions.            Usage:            getAnywhere(x)      argsAnywhere(x)            Arguments:            x: a character string or name.            [...]            > getAnywhere("predict.lda")            no object named 'predict.lda' was found            > library(MASS)            > getAnywhere("predict.lda")            A single object matching 'predict.lda' was found      It was found in the following places      registered S3 method for predict from namespace MASS      namespace:MASS      with value            function (object, newdata, prior = object$prior, dimen, method = c("plug-in",      "predictive", "debiased"), ...)      {      if (!inherits(object, "lda"))      stop("object not of class \"lda\"")      [...]            Il est alors possible de lire le code-source de la fonction predict.lda() et de l'adapter :       # fonction obtenue par predict.lda <- getAnywhere("predict.lda")            ##############################################################################            predict.lda <- function (object, newdata, prior = object$prior, dimen, method = c("plug-in", "predictive", "debiased"), ...) {            ##############################################################################            if (!inherits(object, "lda"))      stop("object not of class \"lda\"")      if (!is.null(Terms <- object$terms)) {      Terms <- delete.response(Terms)      if (missing(newdata))      newdata <- model.frame(object)      else {      newdata <- model.frame(Terms, newdata, na.action = na.pass,      xlev = object$xlevels)      if (!is.null(cl <- attr(Terms, "dataClasses")))      .checkMFClasses(cl, newdata)      }      x <- model.matrix(Terms, newdata, contrasts = object$contrasts)      xint <- match("(Intercept)", colnames(x), nomatch = 0L)      if (xint > 0L)      x <- x[, -xint, drop = FALSE]      }      else {      if (missing(newdata)) {      if (!is.null(sub <- object$call$subset))      newdata <- eval.parent(parse(text = paste(deparse(object$call$x,      backtick = TRUE), "[", deparse(sub, backtick = TRUE),      ",]")))      else newdata <- eval.parent(object$call$x)      if (!is.null(nas <- object$call$na.action))      newdata <- eval(call(nas, newdata))      }      if (is.null(dim(newdata)))      dim(newdata) <- c(1L, length(newdata))      x <- as.matrix(newdata)      }      if (ncol(x) != ncol(object$means))      stop("wrong number of variables")      if (length(colnames(x)) > 0L && any(colnames(x) != dimnames(object$means)[[2L]]))      warning("variable names in 'newdata' do not match those in 'object'")      ng <- length(object$prior)      if (!missing(prior)) {      if (any(prior < 0) || round(sum(prior), 5) != 1)      stop("invalid 'prior'")      if (length(prior) != ng)      stop("'prior' is of incorrect length")      }      means <- colSums(prior * object$means)      scaling <- object$scaling      x <- scale(x, center = means, scale = FALSE) %*% scaling      dm <- scale(object$means, center = means, scale = FALSE) %*% scaling      method <- match.arg(method)      dimen <- if (missing(dimen))      length(object$svd)      else min(dimen, length(object$svd))      N <- object$N      if (method == "plug-in") {      dm <- dm[, 1L:dimen, drop = FALSE]      dist <- matrix(0.5 * rowSums(dm^2) - log(prior), nrow(x),      length(prior), byrow = TRUE) - x[, 1L:dimen, drop = FALSE] %*% t(dm)      dist <- exp(-(dist - apply(dist, 1L, min, na.rm = TRUE)))      }      else if (method == "debiased") {      dm <- dm[, 1L:dimen, drop = FALSE]      dist <- matrix(0.5 * rowSums(dm^2), nrow(x), ng, byrow = TRUE) -      x[, 1L:dimen, drop = FALSE] %*% t(dm)      dist <- (N - ng - dimen - 1)/(N - ng) * dist - matrix(log(prior) -      dimen/object$counts, nrow(x), ng, byrow = TRUE)      dist <- exp(-(dist - apply(dist, 1L, min, na.rm = TRUE)))      }      else {      dist <- matrix(0, nrow = nrow(x), ncol = ng)      p <- ncol(object$means)      X <- x * sqrt(N/(N - ng))      for (i in 1L:ng) {      nk <- object$counts[i]      dev <- scale(X, center = dm[i, ], scale = FALSE)      dev <- 1 + rowSums(dev^2) * nk/(N * (nk + 1))      dist[, i] <- prior[i] * (nk/(nk + 1))^(p/2) * dev^(-(N -      ng + 1)/2)      }      }      posterior <- dist/drop(dist %*% rep(1, ng))      nm <- names(object$prior)      cl <- factor(nm[max.col(posterior)], levels = object$lev)      dimnames(posterior) <- list(rownames(x), nm)      list(class = cl, posterior = posterior, x = x[, 1L:dimen,      drop = FALSE])            } # fin de fonction predict.lda                  Une fois le code lu, il est possible de produire de façon presque lisible les probabilités désirées :             ##############################################################################            probas.lda <- function(object,type="class",newdata) {            ##############################################################################            # ng est le nombre de groupes      # on le connait via les probabilités a priori            ng <- length(object$prior)            # N est le nombre d'individus            nbi <- object$N            # calcul des moyennes pondérées par les probabilités            moyp <- colSums(object$prior * object$means)            # matrice de passage (colonnes LDi)            ldi <- object$scaling      Terms <- object$terms      Terms <- delete.response(Terms)            # données centrées seulement (et sans la constante)            x <- model.matrix(Terms, newdata, contrasts = object$contrasts)      xint <- match("(Intercept)", colnames(x), nomatch = 0L)      x <- x[, -xint, drop = FALSE] }      x <- scale(x, center = moyp, scale = FALSE) %*% ldi            # changement de dimension            dimen <- length(object$svd)      dm <- scale(object$means, center = moyp, scale = FALSE) %*% ldi            # il faut calculer la matrice des distances            dm <- dm[, 1L:dimen, drop = FALSE]      mdist <- 0.5 * rowSums(dm^2) - log(object$prior)      dist <- matrix( mdist, nrow(x),length(object$prior), byrow = TRUE)      dist <- dist - x[, 1L:dimen, drop = FALSE] %*% t(dm)      dist <- exp(-(dist - apply(dist, 1L, min, na.rm = TRUE)))            # voici enfin les probabilités de classement            posterior <- dist/drop(dist %*% rep(1, ng))      nm <- names(object$prior)      cl <- factor(nm[max.col(posterior)], levels = object$lev)      dimnames(posterior) <- list(rownames(x), nm)            return(invisible(posterior))            } # fin de fonction probas.lda            Mais bien sûr rien ne vaut les explications mathématiques. On pourra par exemple consulter la page 20 du document desbois pour s'en rendre compte (copie locale). D'abord par la pratique de la programmation. Le mieux est de faire un maximum d'exercices, d'écrire du code R, mais aussi de lire le code R des fonctions des packages. On pourra utiliser getS3method() pour avoir le code d'une méthode S3, getAnywhere() pour une fonction dite "non visible". Voici un exemple de code R de telles fonctions : library(lattice)      print( methods(histogram) )      print( histogram )      getAnywhere(histogram)      # not run: print( histogram.factor )      getAnywhere(histogram.factor)            Et le résultat de son exécution.       > library(lattice) ;            > print( methods(histogram) )      [1] histogram.factor* histogram.formula* histogram.numeric*            Non-visible functions are asterisked            > print( histogram )      function (x, data, ...)      UseMethod("histogram")      <environment: namespace:lattice>            > getAnywhere(histogram)      A single object matching 'histogram' was found      It was found in the following places      package:lattice      namespace:lattice      with value            function (x, data, ...)      UseMethod("histogram")      <environment: namespace:lattice>            > print( histogram.factor )      Erreur dans print(histogram.factor) : objet 'histogram.factor' introuvable            > getAnywhere(histogram.factor)      A single object matching 'histogram.factor' was found      It was found in the following places      registered S3 method for histogram from namespace lattice      namespace:lattice      with value            function (x, data = NULL, xlab = deparse(substitute(x)), ...)      {      ocall <- sys.call(sys.parent())      ocall[[1]] <- quote(histogram)      ccall <- match.call()      if (!is.null(ccall$data))      warning("explicit 'data' specification ignored")      ccall$data <- environment()      ccall$xlab <- xlab      ccall$x <- ~x      ccall[[1]] <- quote(lattice::histogram)      ans <- eval.parent(ccall)      ans$call <- ocall      ans      }      <environment: namespace:lattice>            Ensuite, il faut essayer de maitriser les outils de la programmation avancée que sont les expressions régulières, les tables de hachage et les tableaux associatifs et la récursivité. Puis il faut s'entrainer à écrire des classes d'objets lorsque les données s'y prêtent, quitte à profiter de packages évolués comme le package iterators qui permet de boucler sur des structures d'objets, comme le package R.oo ou comme le package proto... Enfin (1) , il faut rentrer dans les détails de R : tout y est fonction (et donc objet), même +, [ et <- : `+`(2,5)      "ma variable" <- 2      `<-`(x,3)      y <- 1:5      `[<-`(y,2,-6)      z <- list(a=5,b="oui")      `[[<-`(z,1,-6)      `$<-`(z,"b","non")                  > `+`(2,5)      [1] 7            > "ma variable" <- 2            > `<-`(x,3)            > y <- 1:5            > `[<-`(y,2,-6)      [1] 1 -6 3 4 5            > z <- list(a=5,b="oui")            > `[[<-`(z,1,-6)      $a      [1] -6            $b      [1] "oui"                  > `$<-`(z,"b","non")      $a      [1] -6            $b      [1] "non"                  Enfin (2) , il faut rentrer dans les détails de R et comprendre la notion de fonction de base comme plot, de fonction interne comme plot.xy et de fonction primitive comme `<-`  : plot      plot.xy      `<-`                  > plot      function (x, y, ...)      UseMethod("plot")      <bytecode: 0x33e5380>      <environment: namespace:graphics>            > plot.xy      function (xy, type, pch = par("pch"), lty = par("lty"), col = par("col"),      bg = NA, cex = 1, lwd = par("lwd"), ...)      .Internal(plot.xy(xy, type, pch, lty, col, bg, cex, lwd, ...))      <bytecode: 0x43ac5c0>      <environment: namespace:graphics>            > `<-`      .Primitive("<-")            Enfin (3) , il faut rentrer dans les détails de R et comprendre la notion d'arbre syntaxique, d'évaluation paresseuse, d'environnement, de pile des appels. tmp1 <- quote(if (x>1) { "oui" } else { "non" } )      tmp2 <- as(tmp1,"list")      tmp3 <- lapply( tmp2, typeof )      tmp4 <- objects(pattern="tmp*")            print( tmp1 )      print( tmp2 )      print( tmp3 )      print( tmp4 )            apropos("lazy")                  > tmp1 <- quote(if (x>1) { "oui" } else { "non" } )            > tmp2 <- as(tmp1,"list")            > tmp3 <- lapply( tmp2, typeof )            > tmp4 <- objects(pattern="tmp*")            > print( tmp1 )      if (x > 1) {      "oui"      } else {      "non"      }            > print( tmp2 )      [[1]]      `if`            [[2]]      x > 1            [[3]]      {      "oui"      }            [[4]]      {      "non"      }                  > print( tmp3 )      [[1]]      [1] "symbol"            [[2]]      [1] "language"            [[3]]      [1] "language"            [[4]]      [1] "language"                  > print( tmp4 )      [1] "tmp1" "tmp2" "tmp3" "tmp4"            Enfin (4) , il faut peut-être aussi penser à la gestion des données et savoir utiliser with(), within(), subset(), merge(), stack(), unstack, reshape(), Enfin (5) , il faut peut-être aussi penser aux statistiques et savoir utiliser glm(), nlm(), loglm()... Enfin (6) , il faut peut-être aussi penser à la production automatique de rapports avec R et savoir utiliser toLatex(), RweaveLatex(), Sweave(), de même qu'il faut sans doute s'intéresser aux liens entre XML, les statistiques et R, se dire qu'il est bon de savoir éxécuter du code R pour effectuer des calculs ou tracer des courbes depuis une page Web comme pour valeurs diagnostiques ou aqt, qu'interfacer des bases de données et R en dynamique est important aussi, comme dans la LEAPdb... Bref, il y a encore des milliers de choses à apprendre encore, vu qu'il y a des milliers de fonctions dans des milliers de packages... combien exactement ?    Enoncés de la séance     Retour à la page principale du cours

 

 

   Retour à la page principale de   (gH)