Valid XHTML     Valid CSS2    

 

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)

  1. 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.

  2. 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")     
              
              
    
  3. 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...

  4. 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.

  5. 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...

  6. 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.

  7. 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).

  8. 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 ?

 

 

retour gH    Retour à la page principale de   (gH)