VincentISOZ R
VincentISOZ R
VincentISOZ R
R you ready?
Vincent ISOZ, Daname KOLANI Sciences.ch
Exercice 148.: Comparaison de proportions sur une même population (test binomial
exact) ....................................................................................................................... 592
Exercice 149.: Intervalle de confiance de la proportion .......................................... 593
Exercice 150.: Comparaison de proportions sur 2 échantillons indépendants ........ 594
Exercice 151.: Intervalle de confiance de l'écart-type (test du khi-2 de la variance)
................................................................................................................................. 596
Exercice 152.: Test de Fisher d'égalité des variances ............................................. 597
Exercice 153.: Test de Levene d'égalité de deux variances .................................... 598
Exercice 154.: Robustesse de tests statistiques ....................................................... 600
Exercice 155.: Transformations de Box-Cox .......................................................... 602
Exercice 156.: Transformations de Johnson............................................................ 604
Exercice 157.: ANOVA (Analyse de la Variance).................................................. 606
Exercice 1.: ANOVA à un facteur fixe (ANOVA-1 canonique) empilé ........................ 607
Exercice 2.: ANOVA à un facteur fixe (ANOVA-1 canonique) désempilé ................... 610
Exercice 3.: ANOVA à deux facteurs fixes (ANOVA-2 canonique) sans répétitions avec
ou sans interactions ......................................................................................................... 613
Exercice 4.: ANOVA à deux facteurs fixes (ANOVA-2 canonique) avec répétitions
(réplications) avec ou sans interactions .......................................................................... 622
Exercice 5.: Comparaisons multiples du test de Student avec correction de Bonferroni 627
Exercice 6.: Test de (l'étendue) de Tukey HSD .............................................................. 629
Exercice 7.: Test de Levene et Bartlett d'égalité des variances d'une ANOVA canonique
......................................................................................................................................... 631
Exercice 8.: ANOVA Imbriquée (emboîtée)/Hiérarchique complète............................. 633
Exercice 9.: ANOVA Carré Latin ................................................................................... 636
Exercice 10.: ANCOVA (Analyse de la Covariance) ..................................................... 640
Exercice 11.: MANOVA ................................................................................................ 644
Exercice 158.: ACP (Analyses en Composantes Principales) paramétrique .......... 647
Exercice 159.: Analyse factorielle exploratoire (AFE) ........................................... 653
AFE avec méthode ACP sans rotation! .......................................................................... 654
AFE avec méthode ACP et rotation VariMax! ............................................................... 655
Statistiques non paramétriques .......................................................................... 656
Exercice 160.: M-estimateurs .................................................................................. 656
Exercice 161.: Régression quantile linéaire (globale) ............................................. 657
Exercice 162.: Test d'ajustement du khi-2............................................................... 662
Exercice 163.: Ajustement d'une loi de Poisson par le Khi-2 ................................. 664
Exercice 164.: Test d'indépendance du Khi-deux d'une table de contingence ........ 667
Exercice 165.: Test exact de Fisher ......................................................................... 668
Exercice 166.: Test de McNemar ............................................................................ 670
Exercice 167.: Test d'ajustement du khi-deux avec correction de Yates ................ 671
Exercice 168.: Coefficient de corrélation de concordances des rangs de Kendall .. 673
Exercice 169.: Kappa de Cohen .............................................................................. 674
Exercice 170.: Étude de la cohérence avec l'alpha de Cronbach ............................ 676
Exercice 171.: V de Cramér (mesure d'association)................................................ 677
Exercice 172.: Test de la somme des rangs signés de Wilcoxon pour 2 échantillons
appariés .................................................................................................................... 679
Exercice 173.: Intervalle de confiance de la médiane (via test du signe à 1
échantillon) .............................................................................................................. 681
Exercice 174.: Test du signe binomial (dixit: test de la médiane pour 2 échantillons
appariés)................................................................................................................... 682
Exercice 175.: Test de Mood (test des médianes) ................................................... 684
Exercice 176.: Test de la somme des rangs signés de Wilcoxon pour 1 échantillon
................................................................................................................................. 688
Exercice 177.: Test de la somme des rangs de (Wilcoxon)Mann-Withney pour deux
échantillons indépendants ........................................................................................ 689
Exercice 178.: Test de la somme des rangs signés de Wilcoxon pour 2 échantillons
appariés .................................................................................................................... 690
Exercice 179.: Test d'ajustement de Kolmogorov-Smirnov.................................... 691
Cas à 1 échantillon .......................................................................................................... 691
Cas à 2 échantillons ........................................................................................................ 692
Exercice 180.: Meilleur ajustement au sens de la maximisation de la log-
vraisemblance .......................................................................................................... 694
Exercice 181.: Test de Mantel-Haenszel-Cochran .................................................. 695
Exercice 182.: Test de Grubbs (test des valeurs aberrantes de Grubbs) ................. 697
Exercice 183.: Test de Dixon (test des valeurs aberrantes de Dixon) ..................... 698
Exercice 184.: Test de Friedman (ANOVA canonique non paramétrique par les
rangs) ....................................................................................................................... 699
Exercice 185.: Test de Kruskal-Wallis (ANOVA canonique à 1 facteur non
paramétrique) ........................................................................................................... 700
Exercice 186.: Copulas Gaussiens et Student bivariés ............................................ 702
Régressions ........................................................................................................ 708
Exercice 187.: Régression linéaire univariée par moindres carrés ordinaires (modèle
Gaussien) ................................................................................................................. 708
Effet de levier (leverage) et distance de Cook ................................................................ 720
Exercice 188.: Régression linéaire univariée par moindres carrés ordinaires avec M-
estimateurs ............................................................................................................... 722
Exercice 189.: ANOVA de la régression linéaire (bivariée) à facteurs fixes ......... 724
Exercice 190.: Régression linéaire univariée par moindres carrés ordinaires (modèle
Gaussien) avec plot 2D/3D du modèle sous-jacent ................................................. 727
Exercice 191.: Régression non linéaire univariée par moindres carrés ordinaires
(modèle Gaussien) ................................................................................................... 731
Exercice 192.: Coefficient de corrélation de Pearson et test associé ...................... 734
Exercice 193.: Coefficient de corrélation de Spearman (Spearman rho) ................ 736
Exercice 194.: Régression linéaire multiple (avec ou sans interactions) ................ 737
Exercice 195.: Facteur d'Inflation de la variance (VIF) .......................................... 742
Exercice 196.: Régression pas à pas ascendante ou descendante ........................... 744
Première méthode (descendante package MASS) .......................................................... 744
Première méthode (ascendante package MASS) ............................................................ 745
Deuxième méthode (descendante sans package) ............................................................ 747
Deuxième méthode (ascendante sans package) .............................................................. 749
Troisième méthode (ascendante sans package) .............................................................. 750
Quatrième méthode (meilleur sous ensemble sur R2 ou Cp de Mallows avec package
leaps) ............................................................................................................................... 751
Quatrième méthode (meilleur sous ensemble sur R2, Cp de Mallows, AIC et BIC avec
package biology) ............................................................................................................. 753
Exercice 219.: Carte de contrôle par mesures X barre-S simple ............................. 874
Courbe d'efficacité .......................................................................................................... 876
Exercice 220.: Carte de contrôle par mesures S barre-S simple ............................. 877
Exercice 221.: Carte de contrôle par mesures R barre-R simple ............................ 879
Exercice 222.: Carte de contrôle EWMA avec échantillons ................................... 881
Exercice 223.: Carte de contrôle CUSUM .............................................................. 883
Exercice 224.: Carte de contrôle G.......................................................................... 885
Exercice 225.: Carte de contrôle multivariée T2 de Hotelling ................................ 889
Exercice 226.: Étude R&R (reproductibilité/répétabilité) pour données continues 892
Exercice 227.: Coefficient de corrélation intra-classe ICC1 ................................... 895
Exercice 228.: Étude de l'instrumentation de type I ................................................ 898
Exercice 229.: Plan d'échantillonage par attributs................................................... 900
Méthodes brutes d'inférences ............................................................................ 903
Exercices: Monte-Carlo ........................................................................................... 903
Méthode classique (variables pseudo-aléatoires)............................................................ 903
Méthode de séquence quasi-aléatoire de Sobol .............................................................. 906
Méthode de séquence quasi-aléatoire Latin Hypercube ................................................. 907
Méthode d'acceptation/rejet ............................................................................................ 908
Exercices 230.: Jacknife .......................................................................................... 911
Exercices 231.: Bootstrapping ................................................................................. 913
Décisionnel ........................................................................................................ 915
Exercice 232.: Analyse Procédurale Hiérarchique .................................................. 915
Data Mining (Big Data/Business Intelligence Avancée) .................................. 917
Exercice 233.: Partitionnement ............................................................................... 919
Exercice 234.: OneR (règle unique de classification) ............................................. 920
Exercice 235.: Clustering ID3 (ID-3) ...................................................................... 922
Exercice 236.: Classification Ascendante Hiérarchique (CAH/Dendrogramme) ... 926
CAH avec bagging (bootstrap aggregating) ................................................................... 929
Exercice 237.: K-means (k moyennes en groupe) ................................................... 931
Exercice 238.: K plus proches voisins (k-NN: k nearest neighbours) .................... 937
Exercice 239.: Classification bayésienne naïve ...................................................... 939
Exercice 240.: Classification CART (arbres de régression de classification)......... 941
Exercice 241.: Random Forests (CART Boostrap) ................................................. 950
Exercice 242.: Chi-squared Automated Interaction Detection (CHAID) ............... 953
Exercice 243.: Analyse d'affinité (affinity analysis) ............................................... 961
Exercice 244.: Positionnement multidimensionnel ................................................. 967
Exercice 245.: Analyse discriminante linéaire ........................................................ 971
Exercice 246.: Réseaux de neurones ....................................................................... 976
Exercice 247.: Traiter de gros jeux de données....................................................... 982
Exercice 248.: Hadoop (Cloudera)/R Studio Server ............................................... 987
Séries temporelles (séries chronologiques) ....................................................... 988
Exercice 249.: Génération d'une série temporelle à partir de données brutes ......... 988
Exercice 250.: Extraire des sous-ensembles de séries temporelles ......................... 989
Exercice 251.: Lire des données temporelles .......................................................... 990
Exercice 252.: Décomposition d'une série temporelle par modèle additif et
multiplicatif ............................................................................................................. 993
Introduction
R est un projet GNU libre de traitement des données et d'analyse statistiques (sans AUCUNE
GARANTIE COMMERCIALE quant à la fiabilité des algorithmes) mettant en œuvre le
langage de programmation S créé à l'époque par John Chambers et ses collègues lorsqu'il
était aux laboratoires Bell en 1976. Depuis plusieurs années, deux nouvelles versions
apparaissent au printemps et à l'automne. Il dispose de nombreuses fonctions graphiques.
Le logiciel R est considéré par ses créateurs comme étant une exécution de S, avec la
sémantique dérivée du langage Scheme. C'est un logiciel libre distribué selon les termes de la
licence GNU GPL et est disponible sous GNU/Linux, FreeBSD, NetBSD, OpenBSD, Mac OS
X et Windows (ce qui est un énorme avantage dans l'enseignement et le pratique des
universités et instituts de rechercche). Il représente aujourd'hui l'un des objectifs techniques
majeurs de la communauté GNU.
Il convient de signaler que le 6 Avril 2015 Microsoft a acheté la version entreprise de R
(Revolution Analytics) dont la rapidité d'exécution serait de 20 fois supérieure à la version
standard de R et 42 fois plus rapide que SAS1. Ce qui présage que R avec ses 2 millions
d'utilisateurs va encore plus dominer le domaine de l'analyse statistique de données.
Les étudiants provenant d'universités ont la plupart fait leurs premières armes avec
R et sont donc déjà habitués à cet environnement.
La communauté d'utilisateurs (et donc d'intervenants) semble plus vaste que tous
les autres solutions alternatives d'où une plus grande réactivité sur les forums
(même si le nom du logiciel pose parfois problème avec les moteurs de recherche
puisqu'il ne fait qu'une seule lettre...)
Il contient des packages de haut niveau et souvent dès que de nouveaux modèles
théoriques de statistiques (analyses de données en général) sont publiés, ils sont en
premier disponible sur R pour la communauté académique.
1
http://blogs.technet.com/b/machinelearning/archive/2015/04/06/microsoft-closes-acquisition-of-revolution-
analytics.aspx
2
Si vous étouffez avec SAS... prenez une bouffée d'R !
Dans tout ce livre nous avons fait le choix de travailler avec des fichiers *.csv. Pourquoi ce
choix?:
1. Parce que les fichiers de tableurs comme Microsoft Excel ou autres ne sont pas
pérennes sur les très longues durées (au-delà du demi-siècle)
2. Parce que les fichiers XML sont difficile à produire par les pratiquants dans les
laboratoires qui n'ont pas de formation de programmeur et ne connaissant pas les
normes internationales d'échanges de données.
Le lecteur remarquera que le début du présent document rassemble uniquement des points qui
sont aussi faisables avec MS Excel avec plus ou moins de facilité (parfois plus rapidement,
parfois pas…). Le but est d'abord de montrer les différences entre les deux logiciels comme le
font de nombreux ouvrages avant de passer à des cas très spécifiques à R.
Enfin signalons que l'idée de ce support de cours sera au long terme de traiter des sujets
suivants (et dans l'ordre indiqué) dont les démonstrations mathématiques ont été faites
pendant les cours théoriques:
Un ordinateur est un outil incomparable entre les mains de celui qui sait. Sous les doigts du
Crétin, c'est un revolver manip par un aveugle au milieu de la foule Chester Himes
Remerciements
Cet e-book n'aurait probablement jamais pu voir le jour sans l'apport désintéressé de
nombreux contribuants à ce bijou qu'est R. Je souhaiterais donc exprimer ma plus grande
grattitude et mes plus sincères remerciments à l'équipe qui maintient R et tous les
contributeurs des différentes packages utilisés et cités dans le présent ouvrage. Je ne
manquerai pas d'ajouter leur nom à tous juste ici quand j'aurai un peu plus de temps.
Remarques:
R1. Les premiers exercices sont donc quasiment les mêmes que ceux effectués dans les cours
de statistique théorique, MSP et d'analyse décisionnelle avec MS Excel, Minitab, MATLAB
et SPSS ou lors du cursus d'ingénierie de gestion projets.
R2. Nous ne présentons ici normalement que des outils utilisant des concepts mathématiques
dont j'expose la démonstration mathématique détaillée (démonstrations disponibles dans mon
livre sur les mathématiques appliquées) dans nos cours et qui sont utilisées dans le cadre de
nos activités de consultants. Si cela vient à ne pas être le cas, le titre du sujet est suivi de
l'abréviation WP qui est l'abréviation de: Without Proof.
R3. Nous avons rédigé ce document uniquement pour le fun afin de valider les résultats
obtenus avec MS Excel et surtout pour jouer avec les théorèmes mathématiques présentés
dans nos cours théoriques (démonstrations mathématiques disponibles dans l'ouvrage de
Vincent Isoz sur les Mathématiques appliquées).
R4. Suite à des questions de la part de lecteurs: Non, nous ne sommes pas rétribués par R
pour leur faire de la pub… Des logiciels comme:
XLStat, SPSS, Minitab, SAS, PSPP, Gauss, MATLAB, Statistica, Stata, Medcalc,
StatsDirect, SigmaXL, NumXL, JMP, Weibull++, Design-Expert, PlanExpert, UNISTAT…
font à peu près pareil relativement aux sujets couverts dans ce document. Il fallait juste que
nous fassion un choix... (nous ne pouvons pas passer notre temps à écrire des supports sur
tous les types de logiciels!) et celui-ci s'est porté sur le logiciel utilisé par des écoles dans
lesquelles nous intervenons parfois en tant que consultants. Cependant si des passionnés
veulent reproduire le contenu du présent livre avec leur logiciel de statistique favori qu'ils
n'hésitent pas! Ce serait même fort intéressant de comparer les résultats!! Si nous avons le
temps nous écrirons le même contenu mais avec SPSS (cela étant déjà fait avec Minitab et
partiellement avec Microsoft Excel).
Pour terminer, nous tenons à remercier ici les quelques collègues et clients qui ont bien voulu
nous faire part de leurs remarques pour améliorer le contenu de ce livre électronique. Il est
cependant certain qu'il est encore perfectible sur de nombreux points et qu'il va encore évoluer
puisqu'il y a encore une petite dizaine de sujets dont les démonstrations mathématiques sont
en cours de rédaction pour la prochaine édition du livre de Vincent Isoz les Mathématiques
appliquées et nous montrererons comment les appliquer aussi avec R.
Si vous souhaitez être informé des nouvelles versions majeures de ce document n'hésitez pas
à écrire un mail dans ce sens à Vincent Isoz: [email protected].
PS: Les fichiers utilisés pour les exemples qui vont suivre ne sont fournis qu'aux entreprises,
administrations et écoles faisant appel à nos services de consultants/professeurs.
Courriel: [email protected]
À ce jour interventions dans plus de ~200 entreprises dont 10 du Fortune 500 selon listing
2009 et 3 universités et écoles d'ingénieurs suisses dans des cours de modélisation de bases de
données et simulations stochastiques du risque. Formation de plusieurs dirigeants de
multinationales en one to one.
Accessoirement j'interviens pour des formations sur des logiciels comme MS Project,
MS Visio, MS Access et une vingtaine d'autres dont je délègue l'organisation à des entreprises
spécialisées dans la formation continue en bureautique (niveau licence universitaire et en-
dessous).
Enfin, je conseille aussi vivement à toute personne souhaitant vraiment maîtriser le sujet de
lire mon e-book sur les Mathématiques Appliquées (~4'900 pages).
KOLANI Daname
Domicilié à ce jour à Casablanca (Maroc)
Courriel: [email protected]
J'ai aussi un intérêt pour l'outil informatique notamment les logiciels comme R Project, VBA -
Excel, MatLab, EViews et SPSS…
J'ambitionne me lancer dans le consulting et par là, contribuer à l'optimisation des décisions
financières et à la mise en place de stratégie financière efficiente cadrant avec le modèle
d'affaire des entreprises.
Avertissements
Le contenu du présent support est élaboré par un processus de développement par lequel des
experts en statistiques parviennent à un consensus. Ce processus qui rassemble des
participants bénévoles recherche également les points de vue de personnes intéressées par le
sujet de cet ouvrage. En tant que responsable du présent support, j'assure l'administration du
processus et je fixe les règles qui permettent de promouvoir l'équité dans l'approche d'un
consensus. Je me charge également de rédiger les textes, parfois de les tester/évaluer ou de
vérifier indépendamment l'exactitude/solidité ou l'exhaustivité des informations présentées.
En publiant des textes, il n'est pas dans l'intention principale du présent support de fournir des
services de spécialistes ou autres au nom de toute personne physique ou morale ni pour mon
compte, ni d'effectuer toute tâche devant être accomplie par toute personne physique ou
morale au bénéfice d'un tiers. Toute personne utilisant le présent support devrait s'appuyer sur
son propre jugement indépendant ou, lorsque cela s'avère approprié, faire appel aux conseils
d'un spécialiste compétent afin de déterminer comment exercer une prudence raisonnable en
toute circonstance. Les informations et les normes concernant le sujet couvert par le présent
support peuvent être disponibles auprès d'autres sources que le lecteur pourra souhaiter
consulter en quête de points de vue ou d'informations supplémentaires qui ne seraient pas
couverts par le contenu du présent site Internet.
Et pour rappel:
Normes et validation
Rappelons conformément à ce que nous avons vu dans le cours théorique qu'il est
indispensable pour le chercheur/statisticien/ingénieur professionnel de se baser sur les normes
suivantes (dans l'ordre des plus utilisées) pour son travail et tous les outils dont le présent
support fait l'objet:
ISO 31:2006
Système international d'unités
ISO 3534-1:1999
Vocabulaire et symbole des statistiques
ISO 2602:1980
Interprétation statistique de résultats d'essais - Estimation de la moyenne - Intervalle de
confiance
ISO 3301:1975
Interprétation statistique des données - Comparaison de deux moyennes dans le cas
d'observations appariées
ISO 5479:1997
Interprétation statistique des données - Tests pour les écarts à la distribution normale
ISO 3494:1976
Interprétation statistique des données -- Efficacité des tests portant sur des moyennes et des
variances
ISO 11453:1996
Interprétation statistique des données - Tests et intervalles de confiance portant sur les
proportions
ISO 16269-4:2010
Interprétation statistique des données Détection et traitement des valeurs aberrantes
ISO 16269-6:2005
Interprétation statistique des données - Détermination des intervalles statistiques de tolérance
ISO 16269-8:2004
Interprétation statistique des données - Détermination des intervalles de prédiction
ISO/TR 18532:2009
Lignes directrices pour l'application des méthodes statistiques à la qualité et à la normalisation
industrielle
ISO 3534-3:1999
Plans d'expérience (ou AFNOR NF X 06-080 + NF X 06-081)
ISO 8285:1991
Cartes de contrôle de Shewhart
ISO 10017:2003
Lignes directrices pour les techniques statistiques relatives à l'ISO 9001:2000
ISO 13300:2006
Guide général à l'attention du personnel des laboratoires d'analyse sensorielle
ISO 31010:2009
Techniques d'évaluations des risques
ISO 3951:2006
Règles d'échantillonnage pour les contrôles par mesures
ISO 11095:1996
Étalonnage linéaire utilisant des matériaux de référence
ISO 22514-2:2013
Indices de capabilité
ISO 5725
Précision et fiabilité des mesures en laboratoire (on y retrouve le test C de Cochran et aussi
celui de Dixon)
Bref en gros en 2014 les normes relatives à l'analyse statistique peuvent se résumer avec la
cartographie disponible sur le lien suivant:
https://fr.scribd.com/doc/263063855/AFNOR-Cartographie-Normes-Statistiques
et aussi (Final Guidance for Industry and FDA Staff, January 11, 2002):
Vous pouvez m'envoyer un e-mail pour partager ce que vous avez aimé ou détesté dans le
présent document afin d'en assurer une amélioration continue.
Si vous souhaitez compléter le présent support avec un sujet qui vous tient à coeur et pour
lequel vous avez la démonstration mathématique n'hésitez pas à me contacter. J'intégrerai le
sujet en précisant votre nom et prénom.
Notez que malheureusement, je ne peux pas répondre gratuitement à des questions techniques
d'ingénierie ou de problématique d'entreprise par e-mail pour des raisons professionnelles
évidentes.
E-mail: [email protected]
Liens Internet
http://www.sciences.ch (le site Internet compagnon avec les démonstrations mathématiques)
http://www.r-project.org/
http://www.rstudio.com/
http://cran.univ-lyon1.fr/index.html
http://cran.r-project.org/web/packages/
http://www.itl.nist.gov/div898/handbook/
http://www.cookbook-r.com
http://www.rcommander.com
http://rgraphgallery.blogspot.ch/
Forums
http://www.talkstats.com
http://statistiques.forumpro.fr
http://www.les-mathematiques.net
http://www.developpez.net/forums/f1179/autres-langages/autres-langages/r/
http://forums.cirad.fr/logiciel-R/
Forums
http://freakonometrics.hypotheses.org
http://staffweb.cms.gre.ac.uk/~cd02/EUSPRIG/
Bibliographie
Voici la liste de livres d'une qualité pédagogique et de rigueur extraordinaires que j'ai eu la
chance d'avoir entre les mains et dont je recommande l'acquisition. J'en ai lu beaucoup
d'autres mais qui sont tellement mauvais qu'ils ne valent pas la peine d'être mentionnés:
Le lecteur aura donc compris que je recommande très fortement de compléter la lecture du
présent e-book (non exhaustif concernant R!) par la liste de lecture ci-dessous.
Six Sigma Statistics With Excel and Minitab / 386 pages / Éditions
McGrawHill / Issa BASS
ISBN: 9780071542685
Commentaire: C'est avec ce livre que j'ai été initié à Minitab lors de
ma formation Green Belt chez mon deuxième employeur. Il n'y a
aucun démonstration mathématique et très peu d'équations mais les
exemple donnés ont une orientation industrielle très utile et cela avec
un clarté pédagogique remarquable
ISBN: 053845217X
Installer R
En première année d'Université je considère au 21ème siècle que tout étudiant doit déjà savoir
installer un logiciel seul sur MS Windows ou Mac OS (sinon quoi il peut se faire des soucis
pour sa carrière professionnelle). Donc ici je ne vais pas montrer comment installer R mais
juste présenter le résultat de l'installation excepté pour Linux.
Sous Microsoft Windows 7 (cas le plus courant mais scientifiquement le plus mauvais
choix):
Sous Scientific Linux (RedHat) 6.5 (cas le moins courant mais le meilleur choix car
optimisé par le FermiLab et le CERN pour les applications scientifiques) voici comment
procéder pour installer R. D'abord vous ouvrez une fenêtre de terminal:
et nous allons d'abord nous assurer d'avoir les droits suffisants pour installer R en faisant
quelques manipulations non triviales. D'abord nous ouvrons le fichier de gestion des droits
généraux avec l'éditeur nano:
Et dans l'éditeur il faut rajouter la ligne mise en évidence ci-dessous en rouge (il faut
descendre un peu dans le fichier texte pour trouver l'endroit où mettre cette ligne):
Ensuite, nous écrivons (il s'agit de RStudio Server mais il contient le coeur – core – de la
dernière version de R):
Après 10 à 15 minutes (le temps de télécharger environ 60 modules et de les installer) vous
pouvez vérifier que l'installation soit OK.
Ce qui donne:
Voilà cela donne le squelette de base de R pour Scientific Linux. On doit probablement
pouvoir faire plus simple mais je n'ai pas trouvé (et surtout je n'ai pas le temps de chercher).
Pour la suite nous nous restreindrons au cas le plus fréquent chez les praticiens de l'industrie:
MS Windows (sic!).
Divers R
Changer la langue de l'interface
Par défaut, R s'ouvrira dans la langue que vous avez choisi lors de l'installation. Avec les
versions 3.0.2 et antérieures, il n'existe pas à ce jour de menu permettant de changer la langue
après coup.
Dans les multinationales en Suisse, où l'on travaille souvent en 3 ou 4 langues et il est donc
utile de savoir comment changer la langue de l'interface lorsque l'on rédige une
documentation.
Voici la procédure:
et ensuite vous pourrez le changer avec la commande setwd( ) et vérifier que le nouveau
dossier de travail par défaut a bien été changé:
Nous pouvons aussi voir les fichiers qu'il y dans le dossier de notre choix avec la commande
list.files( ). Par exemple:
avec la commande file.info( ) nous pouvons obtenir certaines informations utiles sur le
fichier:
Mais lorsque vous fermerez et ouvrirez à nouveau R, il aura oublié ce changement. Vous
pouvez vérifier cela en tapant la commande getwd( ):
Dès lors, la méthode consiste à créer un fichier Rprofile dans le dossier de démarrage par
défaut avec la commande suivante au début:
Parfois pour certains tests statistiques quand la p-value est très petite, R renvoie des étoiles (*)
comme nous le verrons lors de notre étude de la régression linéaire univariée (par exemple...):
Nous pouvons lancer l'aide générale ave la commande help.start( ) qui ouvrira la page web
correspondante:
Pour obtenir de l'aide sur une commande qui vous semble mystérieuse, voici comment faire
(du moins quand il y a une documentation...):
Nous pouvons aussi obtenir un rappel (aide) des arguments de certaines commande et
fonctions:
au lieu de:
Le premier dossier est le dossier où R vous avait demandé où installer le premier package (du
moins si vous en avez déjà installé un...) et le deuxième et le dossier des packages par défaut
de R (il y en a une trentaine qui sont déjà installés comme vous pourrez le vérifiant en vous y
rendant).
Pour changer le dossier par défaut d'installation d'abord déterminer votre dossier de travail par
défaut en faisant typiquement la commande getwd( ):
Une fois ceci fait, dans ce dossier créez le fichier .Rprofile avec le dossier par défaut désiré
pour les package:
bingo! C'est très pratique lorsque l'on a plusieurs dossiers avec des packages (personnels, de
notoriété publique ou simplement validé par le CRAN).
Pour obtenir une liste détaillée sur les packages voici une solution simple (nous avons limité
le retour des informations à 5 colonnes mais n'hésitez pas à découvrir par vous même ce qu'il
y au-delà) utilisant la commande installed.packages( ):
et pour mettre à jour rapidement plutôt que de passer par le menu nous utiliserons la
commande update.packages( ):
mais en réalité pour voir tous les data set de tous les packages, vous pouvez utiliser la
commande data(package = .packages(all.available = TRUE)):
Mais alors pourquoi est-ce que je vais majoritairement utiliser mes propres set de données par
la suite vous demandez-vous peut-être???
Simplement parce que mes formations Excel, Minitab, SPSS, Tangra et R ont pour objectif
d'être homogènes... voilà simplement l'explication.
Cependant dans la majorité des cas nous enlèverons les accents et les espaces de la ligne de
titre (légende des colonnes) car cela pose problème à R (puisque non conforme aux normes) et
que cela nous ferait perdre du temps de les renommer avec les commande de R (le temps étant
précieux en formation continue...).
http://datamarket.com/data/list/?q=provider:tsdl
Vous pouvez ensuite choisir un jeu de données en cliquant dessus, en alland dans Export:
et ensuite en cliquant sur Short URL, vous obtiendrez le lien direct d'accès à la source de
données:
et ensuite:
et voilà!
Contrôler la mémoire
Lorsque vous travaillez sur des jeux de données énormes, il va vous falloir maîtriser la
mémoire que R utilise. Nous verrons plus tard quand nous manipulerons des variables
scalaires, vectorielles, matricielles ou autres comment supprimer celles-ci de la mémoire mais
ici nous nous intéressons vraiment à un cas plus général.
R 32 bits est limité 4GB (en pratique c'est plutôt 3GB car l'OS a besoin de
fonctionner)
R 64 bits est limité à 128TB sous LINUX (en pratiques faut trouver de la place pour
mettre tout cela...) et 8TB sous Windows
Pour étendre la mémoire disponible à R ou la limiter, vous allez simplement dans les
propriétés du raccourci de R pour rajouter en Mb la taille désirée:
Ensuite, on peut faire un état des lieux avec la commande du garbage collector qui au passage
détruit les objets inutiles dans la mémoire gc( ):
Sinon on peut aussi changer la mémoire allouée à la volée directement dans R plutôt qu'à
l'ouverture (cependant cela n'est pas toujours permis par le système!):
Mettre à jour R
Mettre à jour R reste à ce jour un peu pénible car il faut aller sur le site web du CRANS,
télécharger, etc. etc. etc.
Il vient ensuite:
Une fois le téléchargement terminé, l'installation classique se fait et une fois celle-ci terminée,
nous avons:
Et ensuite:
Il faut savoir que même si cela installe la nouvelle version, cela ne supprime pas l'ancienne!
Fermer R
Quand vous ferrez du scripting, il vous sera utile de connaître la commande qui ferme R. Soit
vous utilisez alors la commande q( ) seule qui aura l'effet suivant:
soit vous savez d'avance qu'il ne faut pas enregistrer et dès lors vous passez la valeur no en
paramètre:
Fichiers R
Il existe trois types de fichier dans R:
Fichiers *.Rhistory
Nous allons commencer par découvrir le type le plus trivial de fichier qui est le *.Rhistory.
Il s'agit simplement d'un type de fichier qui va sauvegarde l'historique de vos commandes.
Pour voir cela, nous nous ouvrons R et écrivons:
Et il vient:
Si nous fermons R et le rouvrons, nous pouvons recharger les commandes passés en allant
dans Fichier/Charger l'historique des commandes:
Mais rien n'apparaîtra quand vous chargerez votre fichier *.Rhistory. Effectivement pour
retrouver vos commandes, il vous faudra jouer avec les flèches du clavier ou pour
remonter ou descendre dans l'historique des commandes.
On comprend alors mieux l'intérêt de sauvegarder l'historique des commandes dans un fichier
script *.R (voir page 1183 pour rappel).
Fichiers *.Rdata
Les fichiers *.Rdata (fichiers binaires très rapides à charger et à enregistrer) vont stocker
toutes les variables, dates frames, listes, matrices, vecteurs, etc. que vous avez créés ou
modifiés pendant une session de travail R.
Cela peut paraître a priori inutile mais en réalité c'est très utile lorsque l'on fait du calcul
massif (big data).
Effectivement, si vous êtes sur une machine dont la mémoire vive est limitée pour des raisons
de budgets il est alors possible de jouer de façon subtile avec ces fichiers pour libérer
temporairement de la mémoire.
Nous créons une simple variable à laquelle nous affectons une valeur scalaire:
Nous pouvons créer un fichier *.RData en quittant R mais ici le but est d'avoir un contrôle
total sous forme de lignes de commandes (pour intégrer cela plus tard dans des scripts). Pour
enregistrer l'état actuel de la mémoire de R nous écrivons (nous ne spécifions par le dossier
donc il prendra l'espace de travail par défaut actif):
Maintenant nous effaçons cette variable de la session en cours et vérifions qu'elle n'existe
plus:
Si vous utilisez la syntaxe suivante, il va enregistrer dans le fichier *.RData par défaut de
votre installation de R:
Effectivement:
Nous voyons alors vite la possibilité de créer autant de fichiers *.RData que désirés et de les
décharger ou recharger à loisir quand cela est nécessaires (dans les grosses application on aura
typiquement un fichier *.RData pour chaque gros data frame de manipulé).
Petit truc à savoir. Si vous sauvegardez tout un espace de travail et en rechargez tout le
contenu de suite après, cela peut faire gagner parfois un peu plus de 50% de la mémoire
utilisée.
Le but final de R est sur le plus ou moins long terme, lorsque l'on est utilisateur, d'écrire des
scripts qui vont automatiser un travail répétitif de manipulation de données, de calculs et
tests/analyses statistiques (comme le VBA dans la suite Microsoft Office par exemple).
Dès lors, au même titre qu'un utilisateur Excel aura des difficultés à faire du code s'il ne
connaît les structures et opérateur algébrique de base, il en sera aussi souvent de même pour
l'utilisateur R.
Donc nous allons commencer par les mathématiques élémentaires ceci dans le but de préparer
bien plus tard un terrain propice au scripting de haut vol.
Lorsque vous allez faire du scripting, vous aurez parfois besoin de faire des tests logiques
utilisant les valeurs booléennes. Il est donc important de revoir les cas les plus courants.
Nous avons d'abord la syntaxe suivante pour les relations d'ordre les plus courants:
Signalons également le piège suivant à cause de la manière dont fonctionnent les ordinateurs
avec le binaire lorsque nous comparons des calculs avec un résultat étalon:
Pas grand chose à dire ici tellement c'est élémentaire (mais attention! si possible respectez
dans la pratique industrielle la nomenclature de Leszinsky-Reddick pour le nommage de vos
variables):
Autre manière d'écrire cependant parfois plus claire pour les petites classes:
Là encore une fois le but est de faire un petit parallèle avec un tableur. Voyons que nous
pouvons effectivement arrondir les résultats de calculs avec les commandes natives floor( ),
ceiling( ), round( ), signif( ) et trunc( ):
Sinon pour arrondir au 5 centimes près, il faudra passer par la même astuce qu'on utilisait
dans les très vieux tableurs:
Quelques classiques pour afficher les nombres comme le proposent les calculatrices en
utilisant la commande formatC( ):
On retrouvera dans la pratique les nombres complexes dans certains modèles de fractales ou
également dans la transformée de Fourier. Donc vérifions que les éléments de base sont bien
là comme dans n'importe quel tableur:
Voici quelques éléments de la théorie des ensembles qui rappelleront des souvenirs de la
petite école à de nombreux lecteurs:
Comme intersect ne prend que deux éléments en entrée, une manière pour généraliser est
d'utiliser la commande reduce( ):
Il peut arriver dans R que nous ayons parfois besoin de faire des calculs arithmétiques
élémentaires. Voyons un exemple. Ouvrez le fichier CalculsSimple.csv:
Remarquez bien qu'il n'y a pas d'espaces dans les titres et que le séparateur est bien la virgule
et non le ";".
Nous allons d'abord le charger dans R avec la commande read.csv et afficher son contenu:
Maintenant, nus souhaiterions simplement une colonne qui contient la différence des deux.
Pour se faire, nous allons utiliser la commande data.frame( ) qui permet de recréer des arrays
dans la mémoire:
ou pour faire comme dans le livre Minitab avec la valeur absolue en utilisation fonction
abs( ):
Ce qui donne:
Pour le coefficient de variation on peut aussi utiliser la fonction cvar( ) du package BioStatR
au besoin.
Les quantiles étant importants dans la pratique, investiguons un peu quelques syntaxes
possibles:
ou encore pour avoir la catégorie modale textuelle (il existe de nombreuses techniques) en
d'autres termes… une table des fréquences:
Il y a aussi la commande describe( ) du package psych qui permet de se faire une idée du
contenu d'un dataframe:
Il y a aussi la commande native brkdn( ) qui signifie "breakdown" qui donne quelques
informations triviales:
Sinon, toujours en restant dans le domaine des packages, nous avons en utilisant le package
pastecs:
Si le skewness est nul, la distribution est symétrique (cela ne signifiant pas que la
symétrie a lieu sur un pic de la distribution car dans le cas bimodal l'axe de symétrie
peut-être entre les deux valeurs modales symétriques). Si le skewness est positif, la
distribution (la valeur modale/médiane) est penchée à droite (ou il y a des valeurs
extrêmes à droite). Si le skewness est négatif, la distribution (la valeur
modale/médiane) est penchée à gauche (ou il y a des valeurs extrêmes à gauche)
Si le kurtosis est nul (platikurtique) alors l'aplatissement est similaire à celui d'une
distribution Normale. Si la valeur est supérieure à 0 alors la distribution d'intérêt
(leptokurtique) est alors plus haute que celle d'une distribution Normale à moyenne
égale (inversement - mesokurtique - si le kurtosis est négatif bien évidemment)
3
Au moins 11 méthodes de calculer le Skewness. Voir l'étude de Tabor, J. (2010), “Investigating the
Investigative Task: Testing for Skewness - An Investigation of Different Test Statistics and their Power to Detect
Skewness,” Journal of Statistics Education, 18, 1-13. www.amstat.org/publications/jse/v18n2/tabor.pdf
1 n
xi x
3
n i 1
1 3
1 n
xi x
2
n i 1
n
2 1
n 1 n 2
n 1
3/2
3 1
n
et:
1 n
xi x
4
n i 1
1 4
1 n 2
xi x
n i 1
( N 1) ( N 1)1 3( N 1)
2 3
( N 2)( N 3)
N ( N 1) 3( N 1) 2
3 1
( N 1)( N 2)( N 3) ( N 2)( N 3)
4 1 3
2
1
5 4 3 1 3
N
MATLAB -0.647(1),-0.77(2)
Mathematica -0.647(1)
Minitab, SPSS, Excel -0.77 (3)
Commençons par deux grands classiques en utilisant les données suivantes utilisées dans le
cours théorique:
Nous pouvons ensuite représenter l'une ou l'autre des ces deux matrices sous la forme des
ellipses de corrélation de la loi Normale bivariée avec la commande plotcorr( ) du package
ellipse:
Au niveau des légendes ce n'est pas super. Nous pouvons faire mieux:
Voyons les cas classiques que l'on retrouvera énormément de fois même en dehors du
scripting (gestion de data frame, génération de graphiques, etc.). Voici déjà une première série
de commandes utiles utilisant les commandes c( ), sum( ), cumsum( ), diff( ), max( ), min( ),
sort( ), length( ):
Toujours dans les vecteurs regardons que la moyenne arithmétique mean( ) et sa version
élaguée soient bien présentes:
Nous retrouvons bien la même valeur que dans le cours théorique et que dans SAS.
Nous retrouvons bien la même chose que dans le cours théorique mais à la différence de SAS
nous avons du indiquer 0.4 au lieu de 0.2...
Pour le M-estimateur biweight de Tukey nous avons plusieurs packages à choix mais nous
allons nous reporter sur la fonction TukeyBiweight( ) du packaage DescTools:
Nous n'obtenons donc pas la même chose qu'avec SPSS qui donne 8.75 (la joie de l'absence
de normes... a priori).
Mais qui là aussi ne donne pas la même chose que SPSS qui donne 9.09...
Encore une fois, lorsque l'on fera du scripting outre manipuler des scalaires, des entités
algébriques (vecteurs ou matrices) nous aurons aussi à manipuler des textes. Voici quelques
alors quelques commandes classiques de manipulations des textes avec les commandes pour
afficher, concaténer, etc.: print( ), paste( ), cat( ), nchar( ), substring( ), strsplit( ),
toupper( ), tolower( ), gsub( ):
Remarque: Pour mettre la première lettre d'un mot en majuscule ou de plusieurs mots il faudra
passer par une regular expression:
Évidemment certaines commandes que nous venons de voir s'appliquent aussi à des vecteurs
(cas très important puisque les dataframes sont composés de vecteurs!):
celui-ci va prendre de la mémoire pour rien car dans le cas présent nous voyons que nous
avons affaire à des facteurs. Dès lors pour optimiser la mémoire de R, nous pouvons lui dire
qu'il s'agit de facteurs et il associera simplement à chaque texte une valeur 1, 2 (ou plus si
nous avons plus de deux facteurs). Pour faire cela, nous utilisons la commande factor( ) tel
que ci-dessous:
Cela donne alors des variables qualitatives non ordonnées. Si nous avons besoin de définir des
variables qualitatives ordonnées, nous utiliserons la syntaxe suivant (attention! si l'argument
levels n'est pas précisé, c'est l'ordre alphabétique qui sera utilisé):
Pour ceux qui feraient du texte mining, voici quelques fonctions utiles après avoir classé tous
les mots dans des composantes de vecteurs utilisant grep( ):
Pour travailler sur les dates nous allons faire une entorse en supposant que nous avons déjà
étudié l'import de fichier *.csv (voir plus loin pour les détails relativement à ce sujet).
Pourquoi? Simplement parce que c'est la manipulation de nettoyage la plus fréquente
(analyses temporelles oblige...) et donc la plus importante puisque la majorité des pays ne
respectent pas la norme ISO 8601. Bon par la même occasion cela vous montrera comment
(dans l'idée) on modifie toute une colonne d'un data frame aussi pour des textes (l'idée étant la
même).
comme vous pouvez le voir avec la commande as.Date( ), les dates sont converties par défaut
au format normalisé! L'argument de as.Date( ) est donc le format de date tel qu'il se trouve à
l'origine dans le data frame pour que R puisse en comprendre la structure.
Comme nous pouvons le voir R prend les paramètres régionaux de la machine pour afficher le
nom des jours. Nous verrons comment gérer cela avec un package un peu plus loin.
Attention pour le moment où viendra le scripting car lorsque l'on boucle sur des vecteurs
contenant des dates celles-ci sont à nouveau par défaut transformées en nombres:
et comme vous pouvez le voir par la même occasion, l'origine de système de date par défaut
de R est 1970-01-01 mais pouvons à l'aide du paramètre origin changer cela comme nous
voulons.
Nous pouvons aussi transformer une colonne de dates en trimestres avec la commande native
quarters( ):
Maintenant voyons des manipulations triviales sur les dates mais qui semblent nécessiter des
packages (je compléterai au fur et à mesure que je rencontrerai d'autres cas chez mes clients):
Il y a d'autres packages qui ont des options similaires pour ajouter des mois à une date.
Sinon pour gérer des dates dans d'autres langues, on installera le package lubridate:
Pour obtenir les dates de chômées de quelques bourses nous avons la commande holiday*( )
du package timeDate de Rmetrics:
Pour travailler sur les heures et minutes nous allons faire encore une fois une entorse en
supposant que nous avons déjà étudié l'import de fichier *.csv (voir plus loin pour les détails
relativement à ce sujet).
et nous allons nous assurer que R comprenne la colonne DateDeCommande comme étant de
type date et heure sinon quoi nous aurons des problèmes avec l'analyse des séries
chronologiques à l'aide de la commande as.POSIXct( ):
Sinon toujours du package timeDate( ) de RMetrics nous pouvons avoir les heures de
différentes villes de référence pour les fuseaux horaires la finance:
Encore une fois... inutile de préciser l'importance des matrices dans l'analyse financière. Il est
donc important aussi de savoir les manipuler pour plus tard faire du scripting.
Voyons d'abord les trois cas classiques de création utilisant les commandes cbind( ) et
matrix( ) et :
Nous retrouvons donc le même déterminant, les mêmes valeurs propres et les mêmes vecteurs
propres que dans les autres logiciels et que de ce qui a été calculé à la main.
Rappel: Une matrice carrée a toujours des valeurs propres. Par contre si elle n'est pas
symétrique elle aura des valeurs propres complexes.
Lorsque l'on fait du scripting, il peut être utile de compiler des statistiques dans une matrice
avec des noms aux colonnes et aux lignes (la matrice sera alors de squelette pour afficher des
données statistiques diverses et variées).
Voyons un exemple:
Renommer des lignes et colonnes de matrices est aussi utile dans la théorie des graphes pour
avoir une sortie écran plus esthétique. Cependant ce n'est pas tout. Considérons la matrice
d'adjacence suivante vue dans le cours théorique de la théorie des graphes:
Ensuite, avec la commande apply( ), nous pouvons connaître le nombre de flèches entrantes
ou sortantes:
Le deuxième paramètre "1" ou "2" indiquant si nous souhaitons le calcul pour respectivement
les lignes ou les colonnes.
Enfin dans le cours théorique dans le cas des copulas Gaussien, nous avons introduit le
concept de décomposition de Cholesky. Vérifions que nous obtenons bien la même chose que
dans le cours théorique:
Nous construisons la matrice stochastique du graphe comme nous l'avons fait à la main dans
le cours théorique:
Nous retrouvons donc bien les valeurs calculées à la main dans le cours théorique.
Les listes sont considérées souvent comme le type de données le plus complexe de R car on
peut y mélanger des vecteurs, des matrices et des textes en une seule entité (très pratique
plutôt que de renvoyer lors de scripts des variables à n'en plus finir à l'utilisateur!). Voyons
cela via un exemple:
Bon voilà un sujet que nous n'avons traité ni dans le livre sur Minitab, ni dans celui sur SPSS
puisque dans ces deux logiciels, la saisie se fait simplement comme dans un tableur.
Dans R, il peut être utile cependant de savoir comment faire en utilisant la commande c( )
pour générer des colonnes (vecteur) et la commande data.frame( ) pour composer une
matrice::
ensuite pour changer les noms des variables on peut utiliser la commande edit( ) p dans une
feuille similaire à celle d'un tableur:
qui affiche exactement le même tableau mais cette fois avec les valeurs qui sont aussi
modifiables directement.
Vous pouvez aussi utiliser fix( ) pour créer un date frame en partant de zéro:
Pas grand chose à dire car c'est assez simple. Il suffit de lire pour comprendre comment
utiliser les options de base des commandes ls( ) et rm( ) et supprimer ou lister les objets en
mémoire:
Nous pouvons aussi voir la taille d'un objet ou de plusieurs dans la mémoire:
Ce qui donne:
Lorsque vous utiliserez des scripts pris d'Internet et que vous en décortiquerez le
fonctionnement il vous faudra parfois déterminer le typage d'une variable. Voyons quelques
exemples utiles à l'aide de la commande class( ):
Algèbre scolaire
Exercice 22.: Plotter (tracer) des fonctions algébriques
R 3.0.2
Bon rien de bien compliqué pour commencer avec les graphes. Nous allons tracer quelques
polynômes avec des petites améliorations graphiques mineures (on complexifiera au fur à
mesure de notre avancement dans le présent document avec de histogrammes, des doubles
axes, etc.) avec les commandes function( ), plot( ), grid( ) et legend( ):
Nous pouvons obtenir la même chose avec la commande curve( ) en améliorant au passage
un peu le résultat (amélioration que l'on aurait pu appliquer aussi directement à l'exemple ci-
dessus):
Ou une fonction exponentielle avec les deux axes en échelle logarithmique (de façon un peu
basique):
Pour avoir que l'axe Y en échelle logarithmique il vous suffit de mettre de changer l'argument
log (avec un petit exemple clin d'oeil scolaire de l'utilité de l'échelle logarithmique):
ou encore:
Évidemment l'intérêt de savoir faire ce type de graphique c'est typiquement en finance lorsque
l'on connaît explicitement le polynôme du deuxième degré de la frontière efficiente, de
pouvoir aussi tracer la security market line dont l'équation de la droite est aussi explicitement
connue!
Avec la package animation, nous pouvons animer un plot dans un fichier *.gif à l'aide de la
commande savegif( ).
Remarque: Faire des animations est très utile surtout pour les statistiques afin de voir la
convergence des distributions ou dans les plots de cartes pour voir des chemins ou des
densités évoluer.
http://www.imagemagick.org
ce qui donne un image animée *.gif enregistrée dans un dossier temporaire de votre disque:
Pour rester dans les bases, voyons la méthode permet de trouver une racine unique à une
équation supposée avoir des racines en utilisant la commande uniroot( ):
Une fois ceci fait, nous pouvons retourner à notre exemple en chargeant la libraire avec la
command library( ) et en cherchant les racines avec uniroot.all( ) et en les plottant avec la
commande points( ):
R a nativement une commande pour dériver. Voici un exemple de dérivation avec l'estimation
de la dérivée en un point:
Ensuite, pour la dérivée numérique nous allons reprendre l'exemple basique fait dans le cours
MATLAB et comparer la différence... D'abord avec un pas relativement standard:
Nous voyons donc apparaître le même problème de précision que dans MATLAB. Si nous
tatonnons un peu nous verrons que le meilleur pas est de l'ordre de 1e-8 ce qui amène à une
erreur relative de l'odre de 1e-7.
et nous aurons le même problème qu'avant si nous prenons le pas comme valant 1e-16.
Pour l'intégration numérique (il existe un package – Ryacas - pour l'intégration symbolique
mais il nécessite une connexion internet constante et de plus le serveur correspondant n'est pas
toujours en ligne).
Pour un double intégrale la syntaxe est un peu subtile dans R. Voici un petit exemple pour la
loi Normale centrée réduite bivariée:
Nous allons voir ici comment résoudre le petit problème d'algèbre linéaire que nous avons fait
dans le cours théorique sur les systèmes linéaires. Il s'agit simplement d'appliquer les
commandes vues plus haut lors de notre découverte des variables matricielles.
Donc nous retrouvons la même concluions que dans le cours théorique: ce système n'a pas de
solutions.
Nous allons voir par la suite comment gérer certaines fonctions de distributions classiques
dans les petites classes et démontrées dans le cours théorique. Cependant, avant cela, il
convient d'abord de voir une commande particulière que l'on écrit set.seed( ).
Considérons donc d'abord la génération suivante d'une variable aléatoire uniforme comprise
entre 0 et 1:
À chaque fois que vous réexécuterez ces deux commandes, vous n'obtendrez pas le même
graphique. Par contre, si vous utilisez set.seed( ) avec un identifiant, alors dans ce cas le
résultat peut être reproduit:
Ainsi, les deux premiers graphs sont identiques par contre le troisième est différent car il
utilise un autre vecteur aléatoire.
Nous allons voir ici encore un cas trivial de manipulation algébrique. Il s'agit de trouver
l'optimisation d'une fonction univariée en utilisant la commande optimize( ) dont la
spécification d'un intervalle de recherche est obligatoire.
Prenons le cas de la parabole que nous avions obtenu lors de notre étude théorique du modèle
de Markowitz (évidemment c'est un cas trivial mais cela marche aussi pour des fonctions
univariées beaucoup plus complexes):
x, y 0.223,0.00226
Utilisons la fonction baleine à bosse que nous avions utilisée dans le cours théorique lors de
notre démonstration de la méthode d'optimisation de Newton-Raphson:
Nous avons alors en utilisant la commande nlm( ) qui signifie newton local minimum:
Ce qui est effectivement conforme à ce que nous avions vu dans le cours théorique comme le
rappelle ce plot en courbes de niveau ci-dessous:
Nous avons parlé dans la formation des très nombreuses techniques empiriques qui existent
pour chercher des optimum et nous avons démontré deux techniques typiques. Nous allons
voir ici ce que propose R nativement avec la fonction optim( ):
Nous allons utiliser la commande nlmb( ) qui cherche un minimum local avec des contraintes
bornées à l'infini en partant arbitrairement du point:
x, y 1.2,0.6
d'abord pour vérifier que nous retrouvons bien les mêmes valeurs que la commande nlm( ):
Après supposons que nous cherchons le maximum dans la zone mise en évidence ci-dessous:
.....
Nous allons donc ici utiliser R pour résoudre le problème de programmation linéaire étudié en
cours suite à la démonstration mathématique de la méthode du simplexe.
Nous avons vu que cela était relativement laborieux à résoudre à la main, très simple à
résoudre dans les tableurs Open Office Calc et Microsoft. Excel. Eh bien voyons maintenant
cela avec R.
et nous retrouvons donc bien tous les résultats des calculs faits à la main et dans les tableurs
dans le cours théoriques.
Par défaut le package linprog n'accepte que les résultats positifs (raison pour laquelle il n'y
pas d'options pour définir le signe des variables).
Dans ce chapitre, nous allons nous concentrer sur des méthodes se focalisant sur l'utilisation
de l'interface native de R comme à l'habitude.
Ici nous allons vouloir d'abord importer au propre le fichier *.csv suivant qui contient une
petite complication assez courante dans les entreprises:
Ou plus explicitement:
Une fois des opérations/calculs/manipulations effectuées sur le data frame nous pouvons
l'exporter à nouveau en *.csv avec la commande write.table( ):
Nous pouvons aussi associer un commentaire à un data frame pour y mettre une petite
remarque:
Le but ici est simplement de découvrir la fonction file.choose( ) qui injectée dans read.csv( )
dans des plus gros scripts (que ce soit en pur R ou en C# avec R.Net) permet à l'utilisateur de
choisir son fichier avec une boîte de dialogue:
Quand nous commencerons à faire des scripts R complexes pour les domaines des SIG
(systèmes d'information géographiques) ou la finance quantitative, nous nous retrouverons
avec des listes importées d'Internet.
Il est parfois plus prudent de les sauvegarder en local car parfois les serveurs les fournissant
sont en maintenant ou disparaissent tout simplement.
Voyons un exemple d'export en *.csv d'une liste fournie par Yahoo Finance avec le package
quantmod et la commande d'export général à R pour les listes qui est dput( ):
Ce qui donne un fichier *.csv avec un extrait de contenu tel que visible ci-dessous:
Un cas très important dans le domaine financier est de pouvoir aller chercher des données en
temps réel sur des sociétés fournisse ce type de service. La technique n'est pas bien différente
d'avec la technique locale:
Souvent les étudiants demandent comment ils peuvent rapidement copier/coller un jeu
reproductible de leurs données de travail sur un forum pour bien montrer où se situe leur
problème avec le jeu de données réelles.
Pour communiquer rapidement ce petit jeu de données sur un forum, vous pouvez alors
utiliser la commande dput( ):
Le but ici va être de voir une manipulation fréquente dans les laboratoires qui consiste à
fusionner des fichiers découlant typiquement de LabView à l'aide de la commande rbind( ).
Nous allons les fusionner dans R tout en sachant qu'au delà de plusieurs centaines de millions
de lignes R n'arrive plus à gérer la fusion de fichier *.csv (je montrerai une technique pour
cela après).
Vous pourrez alors vérifier dans ce fichier d'export que bien que nous ayons fusionné le
deuxième fichier aussi avec la ligne de titre, celle-ci a été éliminée a priori par R dès la fusion
dans R avec la commande rbind( ).
Si les fichiers sont trop gros (centaines de millions de lignes chacun) il vaut mieux si vous
travaillez sur le système d'exploitation Microsoft Windows utiliser le shell MS DOS avec la
commande suivante qui est extrêmement efficace (je l'ai testée jusqu'à obtenir un fichier final
de 11GB en moins de 3 minutes):
Mais attention au petit piège!!! Il faudra au préalable enlever la ligne de titre à tous les
fichiers excepté le premier (sinon quoi la ligne de titre va se répéter plusieurs fois dans le
fichier final).
Les packages de lecture/écriture XLSX ne marchent pas actuellement avec la version x64 de
R. Il faut exécuter la version x32 et de plus avoir la version x32 du Runtime Java d'installé.
Ensuite le but et d'importer le contenu du fichier suivant xlsxDonnees.xlsx dans un data frame:
PS: Remarquez la dernière colonne n'a pas de titre (légende) pour voir quel va être le
comportement du package dans cette situation.
Nous voyons alors ce que deviennent les espaces initialement dans les titres des colonnes
ainsi que la colonne qui n'avait initialement pas de titre.
Utiliser le presse-papiers
Il est possible de mettre un data frame dans la mémoire (presse-papier) pour pouvoir ensuite
par un simple copier/coller se retrouver avec une un tableau type dans MS Excel (ou autre).
Voici cela avec un exemple d'abord pour des colonnes une par une:
Ensuite, il faut suivant la manière dont vos informaticiens ont installé MS Windows dans
votre organisation aller télécharger la bonne version des drivers ODBC Access de Microsoft:
http://www.microsoft.com/en-us/download/details.aspx?id=23734
ce qui donnait au jour où j'écris ces lignes, la page Internet suivante avec un exécutable à
installer:
Ensuite, nous lançons le gestionnaire ODBC de Windows en 32 bits car le package ODBC ne
marche pas correctement en x64:
Nous prenons ensuite le driver qui correspond au premier fichier que nous voulons utiliser
(nous commençons par un *.mdb et nous verrons ensuite pour le *.accdb):
Nous cliquons ensuite sur Terminer et dans la boîte de dialogue qui apparaît il faut
absolument donner au nom de la source de données le nom du fichier Access lui-même
auquel vous voulez par la suite vous connecter:
et nous allons chercher le dossier où se trouve le Access qui nous intéresse. Nous validons par
OK une première fois:
et nous validons ensuite sur OK et retournons dans R. Nous chargeons comme à l'habitude la
librairie et avec la commande odbcDataSources( ) nous testons si R voit bien notre driver:
Pour la suite nous sommes obligés d'utiliser (du moins au jour où j'écris ces lignes) la version
x32 de R pour pouvoir poursuivre et nous utilisons la commande odbcConnectAccess( )
comme il s'agit d'un fichier Access 2003 ou antérieur et ensuite la commande sqlFetch( ) pour
affichons quelques colonnes et quelques lignes pour voir que tout est ok:
Attention!!! Si le fichier Access est en lecture seule ou ouvert par un autre utilisateur l'accès
par R ne sera pas possible.
La méthode est presque exactement la même pour des *.accdb à part qu'on prendra au début
un driver *.odbc pour les fichiers *.accdb et non *.mdb et que le code change un peu pour car
nous devrons utiliser la commande odbcConnectAccess2007( ) (même si c'est pour des
fichiers d'Access 2010 ou du 2013):
Un cas relativement courant de jeux de données se trouvant sur Internet est ceux qui sont
retrouvent compressés dans des fichiers *.zip.
Avec le contenu:
et nous avons bien récupéré les données comme le montre la commande ci-dessous:
Pour l'exemple nous supposerons que le lecteur aura téléchargé et installé XAMPP:
Une fois téléchargé et installé, si vous travaillez avec une version Pro de MS Windows, il
vous dura arrêter le serveur IIS pour poursuivre:
Ouvez ensuite votre navigateur préféré et tapez http://localhost. Vous devrez alors avoir:
http://dev.mysql.com/downloads/connector/odbc/
Une fois les drivers ODBC de MySQL téléchargés, ouvrez le gestionnaire ODBC de
MS Windows et cliquez sur Ajouter:
et dans la boîte de dialogue qui vient, sélectionnez le driver mis en surbrillance ci-dessous:
Cliquez sur Terminer et vous aurez alors une boîte de dialogue qu'il vous faura remplir
comme indiqué ci-dessous:
Vous pourrez alors choisir la base de données disponible dans le serveur MySQL:
Ici, nous avons pris empiriquement pour exemple la base de donnée installée par défaut avec
MySQL et qui se nomme mysql. Nous validons par OK et revenons alors à:
Ensuite, nous retournons dans R et pour voir la liste des tables de la base de données, nous
utilisons la commande sqlTables( ):
Nous pouvons aussi afficher les colonnes d'une table donnée de la base. Ainsi, avec la
commande sqlColumns( ) nous demandons à R de nous montrer la structure de la table
columns_priv:
Il existe encore d'autres nombreuses possibilités comme de faire des requêtes en SQL
directement sur les tables de votre choix mais pour cela vous pouvez vous référer au PDF du
package RODBC qui contient des exemples très bien fait.
Nous allons donc supposer que Oracle Express 10g est installé:
Pour la démo nous utiliserons le compte système (dont le nom est system et le mot de passe
demandé par l'installateur d'Oracle Express lors de... l'installation).
Une fois donc constaté le démarrage d'Oracle Express, nous ouvrons l'application ODBC de
MS Windows:
Nous cliquons ensuite sur Test Connection. Le nom de la base ainsi que le compte et la
password d'accès ava alors nous être demandé:
De retour dans R, nous chargeons la libraire RODBC, regardons le nombre de tables dans la
base par défaut associée au compte système (nous affichons seulement les dix premières
tables) et nous faisons une requête sur une des tables:
bingo!
Le SQL est un langage extrêmement puissant, bien pensé et très rapide pour traiter et
manipuler des données dans des bases de données relationnelles. Il existe ainsi un package
dans R permettant d'utiliser le SQL sur un ou plusieurs data frames en même temps afin de
pouvoir utiliser les jointures. Ce package est extrêmement utile dans la pratique du big data et
plus rapide que les commandes classiques de R.
Voyons quelques exemples qui seront basés sur le tables d'Oracle Express qui ont été
préalablement exportées au format *.csv.
Il n'est malheureusement pas possible de faire tout ce qui est possible dans Oracle avec du
pour SQL.
Pour plus d'informations avec des cas pratiques basés sur le même jeu de données, le lecteur
pourra se reporter à mon e-book sur le SQL (282 pages):
www.sciences.ch/tiny?url=18
Par rapport à ce que nous allons voir il faut savoir qu'il existe toujours plusieurs manières
d'arriver au résultat. Le but de cet exercice est surtout d'étudier des commandes que vous
retrouverez dans d'autres contextes plus importants, sur de nombreux forums ou livres. Bref
c'est une culture générale qu'il faut acquérir pour être à l'aise.
Les commandes str( ), dim( ) ou encore nrow( ) et ncolumn( ) appliqués à des data frame
sont bien utiles pour avoir un résumé de ce que l'on a:
Voici un exemple très qui nous sera utile (entre autres) beaucoup plus loin quand nous
traiterons du problème du voyageur de commerce utilisant la commande native row.names( ):
Donc là évidemment on peut jouer avec énormément de combinaisons mais dans l'idée voici
un exemple de base avec filtres multiples avec des conditions ET ainsi que OU en utilisant la
commande which( ):
La commande which( ) permet aussi plus simplement de savoir les numéros des lignes où
pour une colonne donnée, certaines valeurs sont satisfaites:
Indiquons aussi la commande subset( ) qui permet de filtrer aussi et que l'on retrouve dans de
nombreux scripts R de data mining:
Ou encore la commande %in% que l'on retrouve dans certains scripts de Finance:
Toujours au niveau des commandes élémentaires de filtrage, voyons les commandes all( ) et
any( ):
- any( ) renvoie une valeur booléenne pour signaler si au moins une valeur correspond au
critère dans le vecteur choisi.
- all( ) renvoie un valeur booléenne pour signaler si toutes les valeurs sans aucune exception
satisfont un certain critère ou non.
Toujours avec notre data frame habituel, nous souhaitons savoir quelle est le prix total avec
rabais le plus proche de 75'000.
Voyons cela en utilisant les commandes sample( ) et les paramètres nrow et replace (dans le
cas ci-dessous nous interdisons à un même individu d'apparaître deux fois dans le tirage):
La commande sample( ) peut être aussi utilisée pour générer des vecteur de textes (utile pour
les arbres binomiaux par exemple):
D'abord nous faisons un résumé simple des secteurs d'activité du jeu de données avec lequel
nous allons travailler et ce en n'utilisant aucune commande du package sampling:
Nous avons donc les proportions représentatives de chaque strate des Secteurs sur l'ensemble
de la population.
Prenons maintenant un même nombre d'individus de chacune des strates (il ne s'agit donc pas
là d'un échantillonage stratitifié à probabilités proportionnelles mais d'un échantillonage
balancé):
Afin d'accélérer le traitement et l'analyse de data frame dans R il peut être judicieux de savoir
supprimer des colonnes.
decreasing = TRUE
Pour faire deux tris (ou plus) dans un ordre différent, voici la syntaxe:
En utilisant la commande duplicated( ) vous pouvez chercher les doublons se trouvant dans
un data frame dans ou une ou plusieurs colonnes.
Il arrivera souvent dans la pratique que l'utilisateur copie/colle des données provenant de
MS Excel sous la forme suivante:
ou avec encore plus de colonnes. Suivant les outils que nous verrons par la suite, il est
nécessaire de mettre les données les unes sous les autres. Certes, un simple copier/coller suffit
mais R a aussi une commande prévue à cet effet.
Nous pouvons ensuite continuer notre travail. Il nous faut pour cela installer le package gtools
comme indiqué ci-dessous:
Nous validons par OK autant de fois que nécessaire. Ensuite, nous retournons à nos
commandes:
Nous souhaiterions faire la même chose qu'avant mais avec une petite subtilité nécessaire
(disons plutôt "préférée" par certains utilisateurs). Il s'agit d'empiler les données mais avec
une variable qualitative de l'appartenance de groupe. Pour cela, nous reprenons le fichier
suivant:
Une méthode possible est similaire à celle d'avant mais avec un petit ajout:
Ensuite, avec la commande tapply( ) que nous retrouvons à de nombreuses reprises nous
pouvons nous amuser un peu:
Pour cela, ouvrez le fichier nous allons utiliser la commande native unstack( ):
Nous allons voir ici quelques techniques élémentaires de synthèse de données. Nous sommes
bien évidemment très loin de ce qu'il est possible de faire beaucoup plus rapidement et
esthétiquement avec le tableur Microsoft Excel (voir va formation d'introduction aux tableaux
croisés dynamiques sur Video2Brain.com).
Voyons une première manière sympathique de synthétiser des données dans un tableau à
double entrée avec des lignes de séparation. Nous installons et chargeons d'abord le package
Epi:
Voici une deuxième manière de synthétiser des données en une table de contingence avec les
commandes table( ) et addmargins( ) et ce sans packages4:
4
La commande xtable( ) que nous verrons lors de notre étude des graphiques permet d'avoir le même résultat
Remarque: Pour extraires les valeurs dans des vecteurs il est possible de transformer la table
en un data frame en utilisant la commande as.data.frame(nom_de_la_table).
Nous pouvons également synthétiser des données (très utile pour des diagrammes à barre ou
des camemberts), avec la commande aggregate( ):
Voyons maintenance ce qu'il est possible de faire avec le package nommé gmodels:
et en utilisant par dessus la commande générique ftable( ) pour formater cela un peu mieux:
Ce dernier cas ressemble aux tableaux croisés dynamique du tableur Microsoft Excel mais en
beaucoup moins flexible, rapide et esthétique.
Voilà une situation dans laquelle nous nous retrouvons fréquemment lorsque nous travaillons
avec des gens utilisant un tableur. Considérons le fichier suivant:
avec la commande reshape( ) du package reshape2 nous allons transformer cela en données
transversales (la syntaxe de la transformation est peu intuitive):
Voyons un cas courant dans les data frame consistant à décomposer le contenu du colonne
contenant des chaînes de caractères en en prenant qu'une seule partie en appliquant la
commande sapply( ) dont la syntaxe est loi d'être intuitive dans le cas présent (il faut au
préalable convertir le typage de la colonne des textes avec la commande as.character( )):
Bon je ne suis pas un fan de ce sujet car je considère que le travail doit être correctement fait
en amont pour éviter ce types de situations (c'est évidemment une opinion très personnelle).
Alors quand nous importons dans R il faut spécifier ce qui est une valeur "vide":
Sinon certaines fonctions statistiques ont directement une option pour gérer ce genre de
situation. Effectivement:
avec is.na( ) nous pouvons aussi savoir combien il y a des valeurs absentes ou non afin de
savoir à quoi nous attendre...:
Dans certaines situations simples plutôt que de faire appel à l'artillerie lourde du SQL pour
associer deux data frame, une solution simple est parfois possible. Elle consiste à faire un
mappage entre deux tables en utilisant la commande native merge( ) de R.
Nous allons faire un mappage et observer le résultat (la conclusion devrait ne pas nécessiter
de commentaires de notre part):
Je ne suis personnellement pas un fan de l'utilisation des graphiques en statistiques (on passe
trop de temps à bricoler le design plutôt que de réfléchir à la robustesse du modèle ou de
chercher pour de meilleurs modèles existants). Donc si jamais j'améliorerai les graphiques
données ci-dessous en fonction des questions des participants et étudiants qui suivent mes
cours.
Au niveau esthétique vous verrez que sommes très loin de ce qu'un tableur comme
Microsoft Excel 2007 ou ultérieur peut faire (ce qui maîtrisent totalement Microsoft Excel me
comprendront). De plus le comportement de R est simplement catastrophique lorsque
l'utilisateur redimensionne la fenêtre contenant un graphique. Cependant, techniquement, R
est plus puissant et plus flexible (dommage que Microsoft n'ajoute pas de nouveaux
graphiques à Excel depuis ces 20 dernières années...).
Pour commencer à faire des graphiques divers élémentaires nous allons essayer d'utiliser
toujours la même source de données, à savoir:
Dans le but d'effectuer quelques graphiques classiques de ce tableau mais que MS Excel par
exemple ne saurait pas effectuer sans un traitement préalable ou sans passer par un tableau
croisé dynamique.
Bon là il n'y a pas grand chose à dire. C'est à usage principalement scolaire mais sinon sachez
que c'est très utile de savoir dessiner des flèches et de mettre des textes sur de plots quels
qu'ils soient avec la commande arrows( ):
Voyons un cas très souvent utilisé en probabilité des flèches (mais simplifié dans le cas
présent):
Au passage, le paramètre pour mettre les nombres de l'axe des ordonnées à l'horizontale est
las:
et si nous voulons la zone de traçage mais sans les axes (remarquez au passage la méthode
pour faire un retour à la ligne dans une légende quelconque en utilisant le \n):
ou avec des axes dans un couleur plus claire que les axes:
Et nous pouvons ouvrir aussi très rapidement une fenêtre toute vide avec la commande
plot.new( ) (cette commande efface le graphique actif en cours!!!)
Mais il existe une autre méthode beaucoup plus intéressant et que nous retrouvons dans des
cas pratiques relativement complexes qui consistent à pouvoir ouvrir autant de fenêtre
graphique que l'on désire avec un nom et une taille au choix:
Ce qui donnera:
Signalons qu'en installant le package plotrix il est possible à tout hasard... de faire des textes
circulaires avec la commande arctext( ):
En entreprise on cherchera à avoir de multiples graphes dans une même fenêtre de R (voir
plus loin) avec des gauges (speedomètres) pour montrer la performance d'un processus, d'un
département ou d'une chaîne de production. Une idée simple est d'utiliser alors une image de
fond comme-ci après en utilisant la commande rasterImage( ) du package png.
http://gastonsanchez.com/blog/how-to/2013/01/10/Gauge-Chart-in-R.html
Lors de la création de scripts vous créerez probablement vos propres graphiques. Dès lors,
dans certains cas il est bien utile de savoir dessiner des rectangles.
Il faudra faire appel au package plotrix pour dessiner des cercles ou des disques avec la
commande draw.circle( ):
Diagramme de Venn
En bioinformatique nous avons régulièrement des vecteurs dont les composantes sont des
séquencages. L'idée est de voir par exemple entre trois bactéries quels sont les dénominateurs
communs au niveau des gènes. Nous pouvons alors typiquement représenter cela avec un
diagramme de Venn. Comme il s'agit le plus souvent d'une représentation de cercles, j'ai jugé
bon de le mettre ici. Petit couac cependant... il doit y avoir une méthode plus simple que ce
qui est proposé ci-dessous... mais à chercher (n'hésitez pas à me dire si vous avez une astuce).
Nous allons devoir d'abord utiliser le package Bioconductor qui n'est pas dans CRAN et une
sous libraire qui se nomme limma:
Laborieux...
Évidemment... Nous ne pouvons pas ne pas montrer comme dessiner des ellipses (pensez ne
serait qu'à cause des fonctions bivariées Normales ou des copulas gaussiens).
Pour cela, nous allons aussi devoir utiliser la package plotrix avec la commande
draw.ellipse( ):
Savoir dessiner des polygones est extrêmement important dans R car dans de nombreux cas
pratiques nous remplissons des surfaces se trouvant en-dessous de la courbe d'une fonction
(particulièrement dans les domaines des probabilités ou des statistiques géographiques où les
territoires sont définis par des polygones).
Donc à partir de là on devine assez aisément la suite (ou voir les exemples dans l'aide de la
commande polygon( )). Par exemple, si nous partons de:
Après, et toujours sur la même idée, on peut superposer plusieurs tracés de lignes remplies les
unes par dessus les autres.
Ou toujours dans les cas simples (nous reviendrons sur les distributions de probabilités plus
loin):
Ou encore pour représenter l'évolution d'un intervalle de confiance (dans le cas présent celui
d'une proportion):
Ou en étant plus exacte puisque ci-dessus nous avons utilisé l'approximation par une loi
Normale:
Une famille de diagramme parfois by sympathique et facilement lisible par le plus grand
nombre en utilisant la commande stripchart( ):
On peut ajouter un bruit pour voir la grossièrement la densité avec l'attribut jitter du
paramètre method:
Commençons donc par le plus simple des graphiques à points plot( ) imaginable avec une
petite couleur de fond différente du blanc (comme cela on apprend deux petites choses d'un
coup):
Avec la commande polygon( ) associée à locator( ) nous pouvons laisser l'utilisateur cliquer
de façon interactive dans le graphique pour créer un polygone:
Ou encore avec la commande native locator( ) seule en cliquant ensuite à différents endroits
du graphe:
Nous obtenons alors les coordonnées des points sur lesquels nous avons cliqués:
Ou encore dans les classiques scolaires de la petite école (attention! pour ploter toute fonction
algébrique dans le plan reportez vous à l'exemple de la page 145):
On peut s'amuser à les relier par des droites (puisque l'index correspond dans le cas présent
aussi à l'ordre des commandes) en utilisant une technique combinant plot( ) et lines( ):
Ce qui donne:
Ou pour revenir à un cas plus simple en ajoutant un ligne de régression linéaire avec la
commande abline( ):
ou encore:
qui donne:
Avec la commande plot( ), nous pouvons aussi avoir une représentation sympathique dont
données Iris que nous avons utilisé dans le cours théorique pour décortiquer les
développements mathématiques de l'analyse en composantes principales et en associant en
plus des couleurs à chaque catégorie!:
Ce qui donne:
Ce qui donne:
Donnera:
Ce qui donne:
Ce qui donne:
Libre à vous ensuite de relier les points par des lignes comme nous l'avons déjà vu plus haut.
Nous allons voir ce type de diagramme juste histoire de le faire car dans la pratique il est très
peu utilisé. Pour cela, nous allons utiliser les données utilisées lors de notre étude théorique de
l'analyse factorielles (mêmes données que celles utilisées dans le support de cours Minitab),
soit:
Bon ici il s'agit simplement de voir que l'on peut obtenir le même type de diagrammes à
points que celui que nous avions fait dans le cours Minitab (mtb) en utilisant le package
plotrix et la commande dotplot.mtb( ):
Donc nous voyons que comme dans Minitab... l'esthétique n'est pas des meilleures...
Sinon un autre grand classique dans les dot plot en se basant sur le fichier suivant:
L'idée du diagramme tournesol est simple: à chaque fois qu'un point se superpose à un autre
de coordonnées identiques, un pétale sous forme de trait est rajouté et le centre de ses
"pétales" successives sont les points superposés eux-mêmes.
Ou encore:
Évidemment comme toujours avec tout diagramme... l'intéprétation peut être délicate (raison
pour laquelle je me rabats toujours si possible sur des indicateurs empiriques et acceptés par
la communauté scientifiques comme étant robustes).
La commande coplot( ) génère en quelque sorte des projections de plot à 3 dimensions sur de
multiples plans de 2 dimensions.
Par exemple, supposons que nous souhaitons comparer les factures payées et non payées en
fonction des articles et de la quantité commandée.
Nous pouvons donc conclure que les factures non payées ne sont pas distinctes pour un article
donné ou une quantité donnée. Les bons ou mauvais payeurs se comportent donc vraiment de
manière aléatoire et non systématique.
Vous pouvez jouer avec ce graphique en remplaçant les factures payées par une variable
numérique ainsi que les noms des articles. C'est très efficace mais la lecture peut vite devenir
compliquée.
Voici la liste des symboles disponibles pour les graphiques à points présentant le type de
symbole pch et leur taille avec cex:
ce qui donne:
Les diagrammes d'associations sont souvent utilisés pour les tables de contingence.
Voici un exemple:
Il faut lire ce graphique dans l'optique d'un tableau de contingence avec la règle du "toute
chose égale par ailleurs" qui nous sera donné par:
Ainsi, toute chose étant égale par ailleurs, les AST Intel 150 non payés sont un peu en-
dessous de la valeur à laquelle on aurait pu s'attendre si toutes les choses avaient été égales
par ailleurs, les AST Intel 200 très en-dessous, les Compaq très au-dessus et les IBM un tout
petit peu au-dessus.
Nous souhaiterions afficher ces résultats sous la forme d'un mosaique de barre. Pour cela,
nous utilisons d'abord la commande xtabs( ) pour faire une synthèse simple et ensuite les
fonctions dotchart( ) et mosaicplot( ) pour les graphiques qui nous intéressent (c'est surtout
ce dernier qui est l'élément d'intérêt!):
Nous allons ici construire un diagramme de type carte à barre représentant le nombre
d'occurrence des secteurs dans tout le tableau. Pour ceci, nous écrivons dans un premier
temps:
Mais les barres d'erreurs manquent!!! Et comme elles ne sont pas natives dans R, voici un
petit script pour palier à ce problème (script tout à fait améliorable):
Nous pouvons également faire comme dans les tableurs des séries (pas facile à deviner)
toujours avec barplot( ):
ce qui donne:
Ce qui donne:
Le lecteur aura remarqué que trier le data frame par jour de la semaine n'influence pas l'ordre
des jours dans le plot (sic!).
Remarquez que nous pouvons obtenir les valeurs détaillées de l'histogramme qui sont implicitement
contenues dans la commande hist( ):
Nous pouvons afficher aussi les quantiles en tant que lignes verticales en ajoutant le script suivant:
1) Elle ne différencie pas les variations positives et négatives par des couleurs
2) Le graphique ne conserve par l'ordre d'origine des données dans les vecteurs il faut donc tricher en
précédent les noms des légendes par des numéros pour avoir l'ordre désiré
3) La fonction ajoute systématique une ligne Total qu'il n'est pas possible de supprimer ni de
renommer
En reprenant le même exemple que dans la formation Microsoft Excel cela donne:
Faisons maintenant un histogramme avec les dates (c'est plus difficile qu'avec de simples
nombres) en utilisant le fichier VentesCleanWithTime.csv qui contient pour rappel:
Un grand classique du monde des tableurs qui se passe de commentaires à part qu'à ce jour le
package ggplot2 avec la commande geom_area( ) ne semblent pas gérer les variables
catégorielles sur l'axe des abscisses (d'où le numéro des jours de la semaine au lieu de leur
nom):
Ce qui donne:
Le but ici va être de vérifier que nous retrouvons la même valeur de l'indice de Gini et la
même courbe de Lorenz que celle calculée à la main dans le cours théorique de Techniques de
Gestion.
et ensuite nous utilisons le package ineq avec les commandes Gini( ) et Lc( ):
Un cas classique afin de voir si nous pouvons faire en gros la même chose que certains
tableurs. Encore une fois c'est une excellente démonstration comme quoi dans R il faut
appliquer des recettes de cuisine car la syntaxe est impossible à deviner intuitivement sans se
référer longtemps à la documentation.
Nous installons ensuite deux packages qui sont lattice et latticeExtra (fallait deviner...):
Avant de créer des graphiques multiples dans une même fenêtre il faut comprendre la
structure de ces dernières.
Mais d'abord, comme les combinaisons de graphiques peuvent prendre de la place, il est bon
de savoir comment contrôler la taille de la fenêtre de plot avant d'y mettre quelque chose:
Voici comment on peut s'entraîner à comprendre les choses en utilisant les commandes par( )
et layout.show( ):
Que l'on peut déjà s'amuser à remplir par quelques plots vides:
Avec une petite complication bien utile pour comprendre les choses:
ou encore (on laisse tomber les marges car tant qu'on y met pas des graphiques, nous ne
verrons pas la différence):
Dans ce qui suivra nous allons retrouver plusieurs fois des combinaisons de graphiques,
voyons un exemple simple:
ce qui donnera:
Maintenant que nous savons faire des combinaisons graphiques (lattices), nous pouvons jouer
avec R pour faire des simili-sparklines (il existe un package pour qui se nomme sparkTable
mais malheureusement à ce jour il ne peut qu'exporter le résultat en un fichier *.html, *.tex ou
*.pdf et non pas l'afficher dans R même).
Voici un exemple que nous laisserons au lecteur le soin de personnaliser ou améliorer (nous
mettons en entrée un vecteur aléatoire si l'utilisateur ne précise rien):
Toujours dans les grands classiques, voyons comment faire de multiples histogrammes de
distributions avec la commande hist( ) et de plus superposés avec transparence:
Nous pouvons aussi nous amuser à combiner un plot avec des histogrammes. Considérons le
jeu de données suivant:
ce qui donne:
Une analyse de la FRE n'a pas le problème de l'analyse par densité où se pose la question des
largeurs des intervalles. Il vaut donc mieux prendre l'habitude de lire les FRE plutôt que les
FDP (fonctions de densité de probabilité).
Nous pouvons cependant faire mieux (je ne suis pas encore arrivé à résoudre le problème du
décalage entre la box plot et la ECDF):
Ce qui donnera:
cela donnera:
Ensuite, faisons un camembert un peu plus intéressant. Résumons d'abord nos données de
ventes par article:
Pour faire un diagramme en anneaux la procédure est très loin d'être intuitive. Il faut dans un
premier temps avoir installé le package ggplot2.
D'abord, nous allons préparer notre jeu des ventes en conséquence pour faire un exemple
sympa:
Ce qui donne:
Ce qui donne:
Pour faire un diagramme d'Ishikawa nous allons d'abord charger le package qcc qui nous sera
très utile dans le chapitre de la maîtrise statistique de la qualité:
Voici un diagramme radial fait pour représenter des quantités sur un intervalle de 24 heures
(le script est inclus dans le framework des participants aux cours).
Ce qui donne:
Le but ici va être de vérifier que nous pouvons reproduire le diagramme de Pareto que nous
avions fait dans R.
Remarque: Pour une analyse ABC il suffit de calculer une colonne en pourcents du total des
fréquences et d'assigner une lettre à ABC à chaque intervalle de notre choix:
Il ne nous faut pas manquer le traditionnel Box Plot avec la commande boxplot( ) et nous y
ajoutons un touche de rug( )...:
Nous pouvons aussi faire varier l'épaisseur des box plot en fonction des effectifs avec le
paramètre varwidth:
ce qui donne:
Les diagrammes en violon essayent de combiner les avantages des boîtes à moustaches et des
estimateurs de la densité locale.
Un grand classique dans les tableurs aussi. Dans R il faut encore une fois utiliser la package
plotrix mais cette fois-ci avec la commande pyramide.plot( ). Voici un exemple repris tel
quel du PDF du package:
Ce qui donne:
Une vieille fonction qui nous est bien connue (sinus cardinal) pour introduire la commande
persp( ) (les plot de surface 3D sont très utiles dans l'industrie où l'on mesure des défauts de
surfaces à quelques dixièmes de microns ou même moins!):
Ou encore en utilisant le package fields pour avoir les dégradés dans le bon sens:
Ce qui donne:
Remarque: Il vaut mieux faire une capture d'écran que d'exporter en tant qu'image car la
qualité de celle-ci est discutable.
Surface 3D animée
Nous retrouvons ici encore une fois le package animation pour exporter une surface d'une
fonction analytique exportée sous forme d'un *.gif animé:
Pour voir cela, nous allons reprendre toujours la fonction sinus cardinal et représenter les
points que nous avons calculés.
Nous allons utiliser le package rgl avec les commandes surface3d( ), spheres3d( ), grid3d( )
et decorated3d( ):
ce qui donne:
Ou encore pour rendre la rotation du diagramme interactive, nous pouvons utiliser le package
rgl avec les commandes plot3d( ), points3d( ), spheres3d( ) et segments3d( ) et enfin
planes3d( ):
Là aussi rien bien de technique. Le but étant de voir si nous arrivons à faire un diagramme à
bulle comme dans n'importe quel tableur (évidemment ce sera encore une fois beaucoup
moins esthétique que le résultat obtenu dans le cours Microsoft Excel).
En reprenant les données du cours Microsoft Excel directement dans la console, nous
écrivons:
Ce qui donne:
Pour ce qui est des plots des surfaces de probabilité cumulée utilisant les méthodes de densité
par noyaux (nous avons vu lors de notre découverte des graphiques comment faire le même
travail pour l'expression exacte de distributions lorsque nous n'utilisions pas la méthode de
densité par noyaux) et les quantiles:
S'il s'agit non pas des quantiles mais du support et que nous ne voulons pas 1 seul intervalle
mais plusieurs, voici la méthode:
Nous allons ici utiliser un package relativement courant qui est la package lattice.
ou:
ou encore:
Nous souhaitons voir ici si nous pouvons reproduire avec R le graphe suivant utilisé dans le
cours théorique de théorie des graphes à partir de sa matrice d'adjacence:
Mais nous avons un message d'erreur. Pas de soucis! L'astuce va consister à définir la source
nous même avec la commande setRepositories( ) et de prendre BioC Software:
et bingo!!!
Nous allons maintenant utiliser la matrice de transition de probabilités définie dans le cours
théorique sur la théorie des graphes:
Une fois ceci fait, nous chargeons les packages plyr, migest et circlize et utilisons les
circos.par( ), chordDiagram( ) et circos.trackPlotRegion( ) comme suite (ce n'est de loin
pas évident à deviner):
La qualtié est mauvais il vaut donc mieux faire un export pour avoir un résultat HD:
Ce qui donne:
Nous faisons la synthèse sur les partenaires car toutes les banques sont listées dans les
partenaires et nous comptons combien de fois chacun apparaît en tant que partenaire ainsi les
plus importantes seront les plus représentatives sur le graphe final:
Nous normalisons cela de façon empirique afin d'avoi un taille convenable sur le graph final
(normalisation faite par essais successifs) et nous construisons à l'aide de la commande
graph.data.frame( ) du package igraph l'objet réseau qui sera utilisé un peu plus loin pour le
graph:
Nous pouvons obtenir la liste des relations avec la commaned E( ) (pour "Edges"):
Nous pouvons faire de même en cercle mais comme il y a des milliers de relations dans le cas
présent ce n'est évidemment guère lisible...:
Nous pouvons ensuite faire quelques petits calculs élémentaires. Par exemple pour chaque
banque, connaître le nombre de liens entrants+sortants avec la commande degree( ):
On peut avoir que les lients entrants (et respectivement sortant en mettant mode="out"):
Nous allons voir ici si avec R nous retrouvons le même dendrogramme que dans le cours
Minitab en utilisant le même jeu de données (rappelez-vous qu'avec Minitab nous avions
obtenu avec la distance euclidienne les mêmes résultats que ceux calculés à la main dans le
cours théorique):
La commande heatmap( ) et le résultat associés sont donnés par (voir page suivante):
Cependant on peut faire mieux que l'horrible commande heatmap( ) qui consiste à utiliser la
commande pheatmap( ) en association avec les packages pheatmap et gplots:
Dans le cours d'ingénierie nous avons étudié en détails la théorie mathématique des plans de
mélange avec mesures aux sommets (plans de mélanges centrés) très utilisés en chimie.
Nous allons voir ici que nous pouvons reproduire graphiquement le mélange qui nous avait
servi d'exemple comme cas pratique de la la théorie.
Et ensuite nous saisissons les données dans la joie en utilisant pour finir la commande
ternaryplot( ):
Dans le cours de mathématique théorique nous avons étudiés la théorie des fractales pendant
un temps significatif en utilisant Maple. Nous allons voir ci-dessous qu'il est tout à fait
possible d'obtenir la même chose avec R.
Ce qui donne:
Ce qui donnera une animation *.gif de 20 images avec comme image finale:
Lorsque vous faites des graphiques vous pouvez certes les exporter manuellement en faisant
un clic droit et choisir une option d'export mais dans la majorité des cas, nous automatiserons
de type de procédure dans des scripts. Voici un exemple:
et si nous ouvrons:
L'export en PDF est superbe! Sachant que l'on peut mettre plusieurs graphiques sur une même
figure de sortie, ainsi que du texte et des résultats de calculs, on imagine aisément comme cela
peut être utile pour automatiser du reporting en entreprise.
Encore une fois, même au niveau de géostatistiques, R ne vaut pas ce que l'on peut faire avec
le tableur MS Excel en termes de rapidités et d'esthétique dans les cas les plus classiques
depuis la version 2010 du tableur en utilisant conjointement Power Pivot et Power Map.
Topographie
Commençons par une trivialité en utilisant la matrice volcano inclue nativement avec R:
Ce qui donnera:
Imaginons maintenant que nous avons une matrice avec le type de terrain indiquant la densité
de verdure. Nous pouvons alors utiliser le script suivant:
Ce qui donne:
Nous allons travailler directement avec le jeu de données intégré à ce même package:
Que l'on peut facilement exporter pour son usage personnel (ceci afin de ne pas resaisir la liste
des pays et régions et sous régions les plus courants) par exemple avec la commande suivante
et en faisant ensuite un coller dans un tableur:
Ce qui donne:
Ce qui donne:
ce qui donne:
Ce qui donne:
Maintenant toujours avec le même package zoomons en perspective sur la région mise en
évidence:
Pas mal... mais nous allons faire un peu mieux en ajoutant un seul paramètre de résolution:
Bon c'est pas mal du tout (là aussi plus besoin de MATLAB pour cela!!!!). Faisons encore
mieux en ajoutant les isoclines:
http://freakonometrics.hypotheses.org/2160
Nous utiliserons le fichier suivant qui contient la liste de 32'250 communes français avec leur
nom et latitude/longitude (le package ne permet pas à ce jour de faire la Suisse):
Et maintenant en prenant que les villes qui comprennent l'expression "-sur-Mer" (comme quoi
les techniques de filtrage que nous avons apprises ne sont pas inutiles!):
http://freakonometrics.hypotheses.org/2308
nous pouvons contrôler les couleurs des régions (ce qui n'est pas faisable à ce jour avec la
Suisse):
Ce qui donne:
Nous allons voir ici quelques possibilité du package ggmap. Commençons par apprendre à
afficher des cartes tout simplement avec différents styles en utilisant les commandes
get_map( ) et ggmap( ):
Nous pouvons également tracer des polygones avec la commande geom_poly( ) et des tracés
avec geom_path( ) ou écrire des textes avec la commande annotate( ) en utilisant par
exemple comme source le fichier suivant:
Nous obtenons:
Nous pouvons indiquer des relevés sur une carte en utilisant à nouveau geom_point( ):
Ou encore afficher des barres de données en des lieux précis avec le fichier de données
suivant:
Ce qui donnera:
On peut aussi s'amuser avec faire des cartes à bulles (remarquez les 0.9 dans le légende qui est
mal compris par R car il s'agit de l'alpha... il faudra le masquer a priori par un rectangle
blanc):
Ce qui donne:
Ce qui donne:
Passons maintenant à un cas très utile dans la meilleure démocratie directe du monde. D'abord
téléchargez sur le site web suivant http://www.gadm.org :
Ensuite, nous utilisons les packages rgdal, maptools, plyr, RColorBrewer, ggplot2 avec la
série de commandes suivantes (merci à Ahmadou Dicko!) en plusieurs étapes (car c'est un peu
long).
1ière étape on charge les packages et on lit le contenu de fichier *.shp5 en utilisant
readOGR( ) et fortify( ):
5
Les fichiers *.shp (shapefiles) sont des fichiers au format ArcView et sont donc de facto le standard pour les
fichiers SIG (bien que ce soit un format propriétaire). La plupart des cartes au monde sont dans ce format.
On génère des données fictives pour l'exemple (on a compris le principe de l'import de *.csv):
Et on trace en cumulant presque tout le savoir cumulé jusqu'ici concernant les cartes:
Ce qui donne:
Nous avons les données suivantes de lieu de domicile et de travail d'individus incluant
différentes colonnes avec les moyens de transport (uniquement les trois premières colonnes
vont nous intéresser):
Nous chargeons les packages nécessaires (déjà connus) et les données d'intérêt (3 premières
colonnes):
Le problème c'est que ce dernier fichier ne contient pas les coordonnées des villes/communes
correspondant aux codes de zones. Pour cela nous allons utiliser le fichier suivant:
Ensuite, nous faisons un mappage classique avec la fonction merge( ) déjà connue ce qui
donne alors:
Ce qui donne visuellement (après une attente certaine dépendant de la puisse de l'ordinateur
car il y a quand même près de 2 millions de chemins/traits à dessiner).
Reste plus qu'à exporter en PDF ou SVG pour avoir une image HD.
Rappelons d'abord comment construire une matrice de distance et montrons ensuite comment
en faire un objet TSP avec la commande TSP( ) et en retirer quelques informations
Pour la suite nous utiliserons un matrice des distances inclue dans le package TSP qui
contient quelques informations de distances euclidiennes entre des grandes villes américaines.
Ce qui donne:
Voyons une petite représentative sympathique que certains apprécieront (à nouveau on peu
faire beaucoup mieux avec un tableur mais bon...).
Considérons les données suivantes d'occupation des salles de conférences d'un petit centre de
centre de formation:
Ensuite, nous allons utiliser le long script de Paul Bleicher comme une boîte noire:
Comme il s'agit d'une manipulation simple à effectuer dans R (afficher les histogrammes de
différentes distributions), nous allons donc commencer par ce premier point.
D'abord signalons les conventions d'écriture de R pour les fonctions de distribution (p pour la
probabilité cumulée, q pour le quantile, d pour la densité et r pour réaliser une valeur aléatoire
au hasard):
Distribution Fonctions
Beta pbeta qbeta dbeta rbeta
Binomial pbinom qbinom dbinom rbinom
Cauchy pcauchy qcauchy dcauchy rcauchy
Chi-Square pchisq qchisq dchisq rchisq
Exponential pexp qexp dexp rexp
F pf qf df rf
Gamma pgamma qgamma dgamma rgamma
Geometric pgeom qgeom dgeom rgeom
Hypergeometric phyper qhyper dhyper rhyper
Logistic plogis qlogis dlogis rlogis
Log Normal plnorm qlnorm dlnorm rlnorm
Negative Binomial pnbinom qnbinom dnbinom rnbinom
Multinomial pmultinom qmultinom dmultinom rmultinom
Normal pnorm qnorm dnorm rnorm
Poisson ppois qpois dpois rpois
Student t pt qt dt rt
Studentized Range ptukey qtukey dtukey rtukey
Uniform punif qunif dunif runif
Weibull pweibull qweibull dweibull rweibull
Wilcoxon Rank Sum Statistic pwilcox qwilcox dwilcox rwilcox
Wilcoxon Signed Rank Statistic psignrank qsignrank dsignrank rsignrank
http://cran.r-project.org/web/views/Distributions.html
Nous allons commencer avec la loi Normale (elle reste la plus utilisée à ce jour) en générant
un vecteur Normal et ensuite en affichant son histogramme et des informations y relatives en
utilisant les commandes rnorm( ) et dnorm( ).
qui donne:
Il est aussi possible de travailler avec des données que typiquement des non statisticiens vont
vouloir ajuster quand même à une loi normale:
Et nous pouvons centrer réduire les valeurs aussi (très utile en finance):
Ce qui donne:
Si nous n'avons que des points dont la distribution est inconnue, nous pouvons utiliser la
commande bkde2D( ) du package KernSmooth pour chercher la meilleure densité 3D
approchée par noyaux locaux:
Ellipse de confiance
Dans un autre style en utilisant le package mnormt et la commande dmnorm( ) associée
nous avons un cas très classique beaucoup utilisé dans le domaine industriel (contrôle
statistique des procédés) dans le cas de côtes bi-dimensionnelles:
Ce qui donne:
Enfin, si nous avons deux variables aléatoires centrées réduites corrélées entre elles, nous
pouvons représenter les ellipses de confiance plus simplement en utilisant les données du
cours théorique:
Ce qui donne:
Nous allons voir à partir de ce point comment représenter graphiquement les diverses
fonctions les plus connues associées à des variables aléatoires d'importance majeure!
Cependant avant de commencer par le premier exemple, il est important je pense de
considérer le cas suivant qui dans un même graphique (certes un peu dense et manquant de
couleurs) rassemble les 4 aspectes possibles d'une variable aléatoire normale (v.a. continue
pour rappel):
Pour la loi uniforme (continue) nous aurons en utilisant les commandes runif( ) et dunif( ):
et là contrairement au support sur Minitab et SPSS nous allons faire les autres distributions
classiques vu dans le cours théorique de niveau BAC puisque les commandes ne sont pas
directement visibles.
Voici notre première v.a. discrète. Là aussi avant de poursuivre considérez le graphique
suivant qui en un seul tenant rassemble les 4 aspect d'une variable aléatoire et remarquez
surtout les arrondis (via l'utilisation de la fonction round( )) pour gérer le fait que le support
concerne les valeurs entières:
Nous avons vu dans le cours théorique que la loi géométrique (loi discrète pour rappel) était
en réalité un simple cas particulier de la loi binomiale négative. Voyons comment tracer aussi
cette dernière en utilisant les commandes rgeom( ) et dgeom( ):
Nous pouvons voir qu'il y a a priori un souci ici car une densité supérieure à l'unité cela n'est
normalement pas possible et le regroupement des barres (intervalles) n'est pas des plus
adéquats...
Il existe encore une fois une autre manière d'avoir un plus joli résultat (mais en réalité un peu
incorrecte puisque cela peut donner l'impression que la fonction de densité est continue) et qui
utilise la commande lines( ) et density( ) cette dernière étant une méthode d'estimation de
fonction de densité par noyaux et par défaut utilisant une gaussienne:
Mais là on voit que la ligne de densité devient vraiment n'importe quoi mais au moins la
densité de l'histogramme est juste! Donc bref il y a un compromis à trouver (sinon on
augmente le nombre de simulation de 125 à 10'000 par exemple...).
Puisque la loi binomiale tend vers une loi de Poisson quand le nombre d'essais devient infini
et que la probabilité devient infiniment petite... il nous faut donc l'étudier!
Pour la loi binomiale (loi discrète pour rappel) je n'ai à ce jour pas trouvé mieux que ce qu'il y
a ci-dessous en utilisant les commandes rbinom( ) et dbinom( ) (c'est un peu particulier car
c'est une loi bien évidemment discrète):
Il existe une autre manière d'avoir un plus joli résultat (mais en réalité un peu incorrecte
puisque cela peut donner l'impression que la fonction de densité est continue) et qui utilise la
commande lines( ) et density( ) cette dernière étant une méthode d'estimation de fonction de
densité par noyaux et par défaut utilisant une gaussienne:
Mais là on voit que la ligne de densité n'est pas vraiment juste non plus mais au moins la
densité de l'histogramme est juste! Donc bref il y a un compromis à trouver (sinon on
augmente le nombre de simulation de 125 à 10'000 par exemple...).
Nous avons vu dans le cours théorique que la loi géométrique était en réalité un simple cas
particulier de la loi binomiale négative (loi discrète). Voyons comment tracer aussi cette
dernière en utilisant les commandes rnbinom( ) et dnbinom( ):
Nous pouvons voir qu'il y a a priori un souci ici car une densité supérieure à l'unité cela n'est
normalement pas possible et le regroupement des barres (intervalles) n'est pas des plus
adéquats...
Il existe encore une fois une autre manière d'avoir un plus joli résultat (mais en réalité un peu
incorrecte puisque cela peut donner l'impression que la fonction de densité est continue) et qui
utilise la commande lines( ) et density( ) cette dernière étant une méthode d'estimation de
fonction de densité par noyaux et par défaut utilisant une gaussienne:
Mais là on voit que la ligne de densité devient vraiment n'importe quoi mais au moins la
densité de l'histogramme est juste! Donc bref il y a un compromis à trouver (sinon on
augmente le nombre de simulation de 125 à 10'000 par exemple...).
Puisque la loi binomiale tend vers une loi de Poisson quand le nombre d'essais devient infini
et que la probabilité devient infiniment petite... il nous faut donc l'étudier!
Pour la loi Poisson (loi discrète pour rappel) je n'ai à ce jour pas trouvé mieux que ce qu'il y a
ci-dessous et utilisant les commandes rpois( ) et dpois( ) (c'est un peu particulier car c'est une
loi bien évidemment discrète):
Là aussi il existe encore une fois une autre manière d'avoir un plus joli résultat (mais en
réalité un peu incorrecte puisque cela peut donner l'impression que la fonction de densité est
continue) et qui utilise la commande lines( ) et density( ) cette dernière étant une méthode
d'estimation de fonction de densité par noyaux et par défaut utilisant une gaussienne:
Mais là on voit que la ligne de densité devient vraiment n'importe quoi mais au moins la
densité de l'histogramme est juste! Donc bref il y a un compromis à trouver (sinon on
augmente le nombre de simulation de 125 à 10'000 par exemple...).
La distribution de Poisson tend vers une loi Normale dont l'écart-type est égale à l'espérance
et que l'espérance (in extenso l'écart-type) deviennent très grand. Mais nous avons déjà étudié
la loi Normale comme tout premier exemple, nous n'allons donc pas y revenir.
Nous avons démontré aussi dans le cours théorique que la loi hypergéométrique était une
généralisation de la loi binomiale (et donc par extension une généralisation de la loi
géométrique). Voyons comment l'obtenir.
Donc pour la loi géométrique (loi discrète pour rappel) en se basant toujours sur la même
approche que les cas précédents mais avec les commandes rhyper( ) et dhyper( ), nous
avons:
Nous pouvons voir encore une fois, qu'il y a encore une fois un souci avec le regroupement
des barres (intervalles) qui n'est pas des plus adéquats....
Il existe (encore une fois...) une autre manière d'avoir un plus joli résultat (mais en réalité un
peu incorrecte puisque cela peut donner l'impression que la fonction de densité est continue)
et qui utilise la commande lines( ) et density( ) cette dernière étant une méthode d'estimation
de fonction de densité par noyaux et par défaut utilisant une gaussienne:
La loi exponentielle (loi continue pour rappel) découle de plusieurs manières de la loi de
Poisson comme nous l'avons vu dans le cours théorique. Voyons donc comme l'obtenir
directement et proprement en utilisant les commandes rexp( ) et dexp( ):
qui donne:
La loi du khi-2 (loi continue pour rappel) découle de la loi Gamma comme nous l'avons vu
dans le cours théorique. Voyons donc comme l'obtenir directement et proprement en utilisant
la commande dchisq( ) et on arrête de faire les histogrammes car on a compris le principe!:
qui donne:
La loi de Student (loi continue pour rappel) découle comme nous l'avons démontré dans le
cours théorique du rapport entre une variable aléatoire normalement distribuée et la racine
carée d'un variable aléatoire distribuée selon une loi du khi-2. Voyons donc comme l'obtenir
directement et proprement en utilisant la commande dt( ):
La loi de Fisher (loi continue pour rappel) découle comme nous l'avons démontré dans le
cours théorique du rapport entre deux variables aléatoires indépendantes distribuée selon une
loi du khi-2. Voyons donc comme l'obtenir directement et proprement en utilisant la
commande df( ):
La loi Log-Normale (loi continue pour rappel) découle comme nous l'avons démontré dans le
cours théorique du logarithme d'une variable aléatoire Normale (cas très courant en finance!).
Voyons donc comme l'obtenir directement et proprement en utilisant la commande dlnorm( ):
ce qui donne:
La loi de Weibull (loi continue pour rappel) serait donc empirique comme nous l'avons
mentionné dans le cours théorique (cas très courant en ingénierie!). Voyons donc comme
l'obtenir directement et proprement en utilisant la commande dweibull( ):
La loi Gamma (loi continue pour rappel) serait donc empirique comme nous l'avons
mentionné dans le cours théorique (cas très courant dans les théories statistiques en général!).
Voyons donc comme l'obtenir directement et proprement en utilisant la commande
dgamma( ):
La loi Beta (loi continue pour rappel) serait donc à peu près empirique comme nous l'avons
mentionné dans le cours théorique (cas très courant dans la gestion de projets et certains
modèles d'analyse bayésienne!). Voyons donc comme l'obtenir directement et proprement en
utilisant la commande dbeta( ):
Voici si jamais une petite commande pour rassembler quelques fonctions de distribution
(ensuite à vous de généraliser en fonction de vos besoins):
Ce qui donne:
Statistiques paramétriques
Nous allons donc voir dans ce chapitre si nous retrouvons les résultats calculés à la main dans
le cours théorique pour voir si les résultats correspondent aux démonstrations mathématiques
détaillées et aussi vérifier s'il y a correspondance avec Minitab, SPSS et MS Excel.
Pour que les choses soient claires, nous ne reviendrons pas dans ce chapitre ni dans celui des
Statistiques non paramétriques sur le problème de la p-value et de ses utilisations fallacieuses
que nous avons longuement développement lors de séances théoriques.
Calculer la puissance d'un test a priori ou a posteriori est important dans le démarche
scientifique. Nous allons ici vérifier que nous retrouvons ou non les mêmes résultats que ceux
calculés à la main dans le cours théorique ou avec Minitab et MS Excel.
Remarque: Le package que nous allons utiliser contient de nombreux calculs de puissance et
de taille d'échantillon, mais au même titre que pour le cours Minitab, nous nous limiterons
seulement à ce que nous avons démontré mathématiquement dans le cours théorique.
pwr.norm.test , n, , H 2
Calculer la taille nécessaire d'un échantillon a priori pour avoir une certaine puissance du test
est aussi important dans le démarche scientifique. Nous allons ici vérifier que nous retrouvons
ou non les mêmes résultats que ceux calculés à la main dans le cours théorique ou avec
Minitab et MS Excel.
Pour cela, nous allons encore une fois utiliser la commande pwr.norm.test( ) du package
pwr:
pwr.norm.test , P, , H 2
Nous retrouvons bien les même résultats qu'avec Minitab (et donc les mêmes conclusions)
simplement affiche automatiquement la courbe de puissance ce qui est bien pratique. Nous
pouvons cependant la construire nous même rapidement avec R:
Bingo!
Calculer la puissance d'un test a priori ou a posteriori est important dans le démarche
scientifique. Nous allons ici vérifier que nous retrouvons ou non les mêmes résultats que ceux
calculés à la main dans le cours théorique ou avec Minitab et MS Excel.
Nous retrouvons donc bien les mêmes valeurs que dans Minitab et presque les mêmes que
celles calculées à la main dans le cours théorique.
N'ayant pas trouvé de fonctions dans quelque package que ce soit pour la résolution du test Z,
nous allons donc passer directement au t et nous sauterons le calcul de la résolution du test p
n'ayant pas trouvé de package intégrant aussi la fonction.
Nous retrouvons bien la différence de 2/3 (à une petit erreur d'arrondi près).
Calculer la taille nécessaire d'un échantillon a priori pour avoir une certaine puissance du test
est aussi important dans le démarche scientifique. Nous allons ici vérifier que nous retrouvons
ou non les mêmes résultats que ceux calculés à la main dans le cours théorique ou avec
Minitab et MS Excel.
Pour cela, nous allons encore une fois utiliser la commande pwr.t.test( ) du package pwr:
Nous retrouvons donc bien les mêmes valeurs que dans Minitab et que celles calculées à la
main dans le cours théorique au dixième près.
Nous souhaitons ici simplement vérifier que nous obtenons les mêmes résultats que ceux
obtenus dans MS Excel avec la méthode Monte Carlo (et in extenso vérifier les résultats
obtenus lors de l'étude théorique du test de d'Anderson-Darling) et Minitab.
En utilisant la commande ad.test( ) et avec les mêmes données brutes que dans le cours
théorique, nous avons en utilisant le package ADGofTest:
Donc ici nous voyons un énorme avantage par rapport à Minitab déjà: nous pouvons choisir la
distribution selon nos désirs!
Nous voyons cependant une différence significative entre R et Minitab. Mais nous en
connaissons l'origine puisque nous avons fortement critiqué certains aspect du test d'AD
pendant le cours théorique (cependant la conclusion reste la même).
Il existe un autre package spécialisé des les tests de Normalité et qui est nortest mais bon
voilà il souffre d'un léger problème comme en atteste la capture d'écran ci-dessous:
Accessoirement, rappelons que nous avons vu dans le cours théorique que le test d'ajustement
de Cramer-Von Mises est un cas particulier du test d'Anderson-Darling (il n'y a simplement
pas de dénominateur). Ce test existe aussi sous la commande cvm.test( ) dans le package
nortest mais il souffre du même défaut. Effectivement:
Nous souhaitons ici vérifier que nous obtenons les mêmes résultats que ceux obtenus dans
MS Excel avec la méthode Monte Carlo (et in extenso vérifier les résultats obtenus lors de
l'étude théorique du test de Ryan-Joiner) et Minitab.
Rappelons que l'avantage de ce test non paramétrique de normalité (donc basé sur les rangs)
est sa simplicité mais qu'il n'est pas contre pas adapté lorsque que trop de valeurs identiques
se répètent.
En utilisant la commande shapiro.test( ) et avec les mêmes données brutes que dans le cours
théorique, nous avons:
La sortie est complément différente de Minitab au niveau esthétique (c'est bien dommage!)
mais la valeur W et la p-value sont conformes en ordre de grandeur (voir mon livre Minitab
pour les valeurs numériques exacte à comparer).
Curieusement le test Z n'est pas intégré par défaut de R. Certes il est extrêmement simple de
récrire la fonction correspondante diront certains (et c'est vrai) mais pour moi ce n'est pas un
argument car dans ce cas on peut aussi prendre le temps de récrire toutes les tests statistiques
car ils sont tous simples! Donc soit on a un outil qui nous fait gagner du temps systématique
pour tout soit pour rien mais pas l'entre deux. Pour cela, par exemple, je préfère de loin
Minitab.
Mais le test Z existe quand même dans le package TeachingDemos mais seulement dans le
cas où les données mesurées sont connues (et non pas juste le résumé).
Voyons cela avec le même exemple que dans le cours Minitab, SPSS et MS Excel en utilisant
la commande z.test( ):
Nous retrouvons donc les valeurs de Minitab mais en plus précis et nous voyons in extenso
que les résultats avec MS Excel étaient peu précis.
Il est connu dans un État que les enfants d'un certain âge ont un poids de 45 kilogrammes et
un écart-type de 13 kilogrammes (espérance et écart-type de la population). Un plainte est
posée par des parents d'élèves comme quoi les enfants d'une école sont sous-alimentés. Pour
cela les parents d'élèves s'appuient sur le fait que 25 enfants du même âge ont un poids moyen
de de 40.5 kilogrammes.
Nous pouvons donc bien constater que nous obtenons les mêmes résultats que ceux calculés à
la main et que dans des logiciels comme Microsoft Excel et Minitab.
Le but de cet exercice est plus pédagogique qu'autre chose. Imaginons effectivement que vous
connaissez les estimateurs de la moyenne et l'écart-type d'une loi Normale et la taille de
l'échantillon que vous allez mesurer.
La question pourrait être alors simplement la suivante: si je tirais au hasard des individus
d'une telle population (simulable par du Monte Carlo en théorie), quelle est la probabilité
cumulée (ou 1 – la p value) que ma mesure dépasse une certaine valeur seuil.
Nous avons alors dans le cas particulier unilatéral la possibilité d'utiliser la fonction
as.randtest( ) du package ade4. Ce qui donne:
Là encore nous allons vérifier la conformité des résultats avec toujours les mêmes logiciels et
les calculs effectués à la main suite à la démonstration mathématique faite dans le cours
théorique.
Le test de Student à 1 échantillon est lui intégré par défaut dans R mais là encore seulement si
les données brutes sont disponibles (encore une fois c'est dommage mais bon...).
et là nous sommes en conformité avec Minitab (mais ce dernier ne donne pas la p-value) et
MS Excel (avec qui on retrouvait toutes les valeurs).
Nous retrouvons donc les bons résultats avec les mêmes conclusions.
Comme d'habitude nous allons contrôler que nous obtenons les mêmes résultats que ceux
calculés à la main ou avec MS Excel et Minitab après la démonstration mathématique du test
dans le cours théorique. Nous utiliserons les mêmes données de mesures et la commande
t.test( ):
Là encore nous allons vérifier la conformité des résultats avec toujours les mêmes logiciels et
les calculs effectués à la main suite à la démonstration mathématique faite dans le cours
théorique.
Les résultats sont donc conformes à la différence que la p-value est extrêmement précise dans
R...
Nous retrouvons donc les bons résultats avec les mêmes conclusions.
Nous continuons donc avec le contrôle de conformité avec les démonstrations mathématiques
faites en cours et les calculs faits à la main et dans Minitab et toujours avec les mêmes
données.
Nous retrouvons donc les mêmes valeurs que Minitab et MS Excel. À la différence que
Minitab donne des informations que R ne donne pas et réciproquement (alors qu'avec
MS Excel on pouvait bien évidemment avoir toutes les informations).
Nous retrouvons donc les bons résultats avec les mêmes conclusions.
Nous continuons donc avec le contrôle de conformité avec les démonstrations mathématiques
faites en cours concernant le ratio de deux moyennes et son intervalle de confiance.
En utilisant les mêmes données que pour l'exemple précédent, nous obtenons:
Ce qui corresond à peu près aux calculs faits à la main où nous avions obtenu:
2
X X a
g g 1 g
Y Y b
1,2 0.987,1.084
1 g
Nous continuons donc avec le contrôle de conformité avec les démonstrations mathématiques
faites en cours et les calculs faits à la main et dans Minitab et toujours avec les mêmes
données.
Nous retrouvons donc à peu près les mêmes valeurs que Minitab et MS Excel. À la différence
que Minitab donne des informations que R ne donne pas et réciproquement (alors qu'avec
MS Excel on pouvait bien évidemment avoir toutes les informations).
Nous retrouvons donc les bons résultats avec les mêmes conclusions.
Calculer la taille nécessaire d'un échantillon a priori pour avoir une certaine puissance du test
est aussi important dans le démarche scientifique. Nous allons ici vérifier que nous retrouvons
ou non les mêmes résultats que ceux calculés à la main dans le cours théorique ou avec
Minitab et MS Excel.
Pour cela, nous allons encore une fois utiliser la commande pwr.p.test( ) du package pwr:
Nous retrouvons donc bien les mêmes valeurs que dans Minitab et toujours la même
différence qu'avec ce que nous avons calculé à la main dans le cours théorique.
Calculer la taille nécessaire d'un échantillon a priori pour avoir une certaine puissance du test
est aussi important dans le démarche scientifique. Nous allons ici vérifier que nous retrouvons
ou non les mêmes résultats que ceux calculés à la main dans le cours théorique ou avec
Minitab et MS Excel.
Imaginons que nous souhaiterions mettre en évidence une différence de 20% dans une étude
(ou sondage) où nous nous attendons à avoir 50% de proportion expérimentale (donc une
proportion de comparaison de 50-20=30%)
Pour cela, nous allons encore une fois utiliser la commande pwr.2p.test( ) du package pwr:
Nous retrouvons donc bien les mêmes valeurs que dans Minitab et toujours la même
différence qu'avec ce que nous avons calculé à la main dans le cours théorique.
Nous allons voir ici si nous retrouvons encore une fois les mêmes valeurs que celles calculées
dans le cours théorique à la main ainsi que dans Minitab.
Nous utilisons la commande native poisson.test( ) et en bilatéral pour déterminer pour une
compagnie d'aviation ayant eu 2 deux crashs en 1'000'000 de vols (événement très rare),
quelle est l'intervalle de confiance en bilatéral à 95% sachant qu'au niveau mondial le nombre
d'accidents par millions est de 0.4.
Donc là nous retrouvons les mêmes valeurs que dans le cours Minitab globalement, mais par
contre pour la borne de gauche de l'intervalle, comme nous l'avons déjà mentionné dans le
cours Minitab, le résultat ne correspond pas avec nos calculs faits à la main. Nous avions
cependant montré dans le cours théorique quels étaient les calculs effectués en réalité dans
l'algorithme (cependant sans en trouver la démonstration).
Maintenant procédons à l'exemple unilatéral fait aussi dans le cours Minitab avec les non-
conformes.
Une société fabrique des télévisions en quantité constante et a mesuré le nombre d'appareils
défectueux produits chaque trimestre pendant les dix dernières années (donc 4 fois 10
mesures = 40 trimestres). La direction décide que le nombre maximum acceptable d'unités
défectueuses est de 20 par trimestre et souhaite déterminer si l'usine satisfait à ces exigences
Nous reproduisons comme à l'habitude l'exemple du cours théorique et fait avec Minitab et à
la main:
Une compagnie d'aviation a eu 2 deux crashs en 1'000'000 de vols (événement très rare). Une
autre compagnie a eu 3 crashs en 1'200'000 vols. Quel est l'intervalle de confiance en bilatéral
à 95% en supposant que la différence devrait être nulle.
R n'implément pas cette méthode à ce jour avec les hypothèses sous-jacente du modèle que
nous avons vu dans le cours théorique (hypothèse de continuité, stabilité de la loi de Poisson).
Nous continuons donc avec le contrôle de conformité avec les démonstrations mathématiques
faites en cours et les calculs faits à la main et dans Minitab et toujours avec les mêmes
données.
Je n'ai pas trouvé de package faisant ce calcul (le dernier qui le faisait ayant disparu ce qui est
bien dommage) donc voici une fonction qui fait le test et qui redonne le même résultat que ce
que nous avons calculé à la main, dans MS Excel et Minitab dans le cours théorique:
Nous retrouvons donc les bons résultats avec les mêmes conclusions.
Nous continuons donc avec le contrôle de conformité avec les démonstrations mathématiques
faites en cours et les calculs faits à la main et dans Minitab et toujours avec les mêmes
données.
Nous retrouvons donc exactement les mêmes données que celles calculées à la main et dans
MS Excel et Minitab (nous parlons alors parfois de test binomial exact de Clopper-Pearson).
Nous retrouvons donc les bons résultats avec les mêmes conclusions.
Nous continuons donc avec le contrôle de conformité avec les démonstrations mathématiques
faites en cours et les calculs faits à la main et dans Minitab et toujours avec les mêmes
données.
Dans les deux cas, nous voyons que nous ne retrouvons pas le même résultat que les calculs
faits à la main et dans MS Excel pour l'intervalle de confiance (mais pour le reste c'est OK...).
Il y a une différence d'environ 10% en ce qui concerne les bornes c'est donc non négligeable.
Nous continuons donc avec le contrôle de conformité avec les démonstrations mathématiques
faites en cours et les calculs faits à la main et dans Minitab et toujours avec les mêmes
données.
Doc nous ne retrouvons évidemment pas la même chose que dans Minitab et dans MS Excel
puisque nous y avions fait une approximation par la loi Normale alors que R fait une
approximation par la loi de Khi-2 (bon la conclusion reste toutefois le même dans ce cas
particulier).
Si nous faisons comme Minitab immédiatement un test exact de Fisher (Minitab le fait
automatiquement), nous avons avec R en utilisant la commande fisher.test( ):
Nous voulons déterminer l'intervalle de confiance de cette variance (c'est comme-si nous
cherchions à comparer avec un ratio où nous avons 1 au dénominateur) en utilisant la fonction
sigma2.test du package sigma2tools.
Un échantillon de 9 vis a été tiré d'une ligne de production et la mesure de leur diamètre en
[mm] a été reportée ci-dessous:
et tout est parfait! Nous retrouvons les mêmes résultats que dans le cours théorique et qu'avec
MS Excel et Minitab.
Nous continuons donc avec le contrôle de conformité avec les démonstrations mathématiques
faites en cours et les calculs faits à la main et dans Minitab et toujours avec les mêmes
données.
Nous retrouvons donc à peu près les mêmes valeurs que Minitab et MS Excel. À la différence
que Minitab donne des informations que R ne donne pas et réciproquement (alors qu'avec
MS Excel on pouvait bien évidemment avoir toutes les informations).
Nous continuons donc avec le contrôle de conformité avec les démonstrations mathématiques
faites en cours et les calculs faits à la main et dans MS Excel (Minitab n'ayant pas le test de
Levene mais de Brown-Forsythe même s'il indique le contraire) et toujours avec les mêmes
données.
Nous utilisons donc la commande leveneTest( ) du package car en préparant les données
avec une structure peut triviale:
Nous retrouvons donc exactement les mêmes valeurs que celles calculées dans MS Excel.
Voyons maintenant avec la variante robuste de Brown-Forsythe pour voir si là encore nous
retrouvons les mêmes résultats (et bien évidemment avec les mêmes conclusions). Nous
allons devoir utiliser la commande levene.test( ) du package lawstat:
Et nous retrouvons donc bien la même chose à la différence que R donne la statistique et la p-
value avec plus de précision (comme à chaque fois quoi!).
Nous allons voir ici comme il est facile de mettre en oeuvre avec R ce que nous avons fait
dans le cours MS Excel et qui consistait à tester la robustesse de certains tests simples ou
particulièrement complexes.
Nous n'allons pas ici refaire un exemple d'application pour chacun des tests (qu'ils soient
paramétriques ou non paramétriques comme nous l'avions fait dans le cours MS Excel) car ce
serait probablement plus ennuyeux qu'autre chose pour le lecteur.
Afin de gagner du temps aussi de mon côté (car je déteste recréer la roue excepté si ce n'est
pour faire mieux), j'ai reprise l'exemple mot pour mot, lettre par lettre du e-book de
Vincent ZOONEKYND dont le but est de tester la robustesse du test de Student à la non
normalité des données (avec une distribution uniforme au lieu de Normale) et pour lequel il
faudrait alors normalement utiliser un test non paramétrique de type Wilcoxon U.
Donc nous y allons pour faire 1'000 tests de Student pour échantillons indépendants en
bilatéral et faire la moyenne des p-value qui sont au-delà de 5% (soit les vrais positifs dans le
cas présent):
Donc nous avons 94.9% de vrais positifs. Donc le test est robuste pour une distribution
uniforme qui est non pathologique relativement à la loi Normale. Ent tout cas cela donne
l'approche de la méthode pour généraliser ensuite à tout autre test.
Nous pouvons faire de même avec l'intervalle de confiance (observez bien l'approche
sémantique car c'est très pertinent!):
Ici nous allons juste contrôler que nous retrouvons le même résultat de transformation
(empirique) de Box-Cox que dans le cours Minitab. Nous utilisons pour cela un fichier
contenant 125 mesures dans une seule colonne et utilisons la commande powerTransform( )
du package car:
Donc cette commande suggère que peut normaliser cette variable en la mettant à la puissance
0.0318. Nous retrouvons donc la même valeur que dans Minitab à la différence que R est
beaucoup plus précis et donne la p-value de l'utilité de la transformation.
Nous pouvons ensuite utiliser notre connaissance des graphiques pour observer le résultat
qualitativement:
Bon le but ici ne sera pas de débattre de la justesse de ces transformations (ou de leur origine
scientifique) qui ont donc pour rappel de "normaliser"6 des données. Ceci a aussi déjà été fait
en cours et de plus, il suffit de Googler un peu pour se faire sa propre opinion sur cette
technique "d'ingénierie statistique".
Les transformations proposées sont les mêmes que pour Minitab mais avec une notation
différentes pour les coefficients. Donc si jamais voici un extrait du package:
Nous reprenons donc les mêmes données que pour la transformation de Box-Cox et utilisons
la commande JonhsonFit( ) du package SuppDists:
6
Dans le sens: faire que la distribution ressemble le plus possible à une loi Normale
Bon ceci étant dit voyons ce que nous pouvons tirer de cet ajustement:
Avant de commencer juste en petit rappel sur la notation dans R des modèles de régression
s'avère après expérience indispensable:
Syntaxe Lecture
Z ~ X Y Le modèle est supposé du type Zi X i Yi
Z ~ X *Y Le modèle est supposé avec interactions du type Zi X i Yi X iYi
Z ~ X /Y Le modèle est supposé avec interactions imbriqué Zi X i X iYi
Z ~ X :Y Le modèle est seulement avec interactions Zi X i Yi
D'autres exemples suivront si besoin…
Imaginons une entreprise faisant les trois huit. Nous avons trois équipes qui travaillent sur une
même machine. Nous souhaitons vérifier avec un seuil de confiance de 95% s'il y a une
différence de productivité moyenne entre les trois équipes sur une semaine de travail.
Nous allons vérifier que nous retrouvons bien les valeurs calculées à la main dans le cours
théorique ainsi que dans le cours MS Excel et Minitab en utilisant les mêmes données:
Déjà avant de sortir l'artillerie lourde une analyse qualitative peut parfois suffire:
Donc une analyse plus poussée se justifie bien! Nous utilisons alors la commande native
aov( ) de R:
La commande plot.design( ) nous permet aussi d'afficher les effets moyens sous une autre
forme que les box-plot:
Pour le cas où les données sont désempilées, considérons l'exemple vu dans le cours théorique
et que nous retrouverons lors de notre étude de la régression:
Nous devons donc restructurer les données pour R avec la fonction stack( ):
Le résultat bien que peu explicite est cependant bien conforme aux calculs effectués dans le
cours théorique.
Imaginons une entreprise faisant les trois huit. Nous avons trois équipes qui travaillent sur une
même machine. Nous souhaitons vérifier avec un seuil de confiance de 95% s'il y a une
différence de productivité moyenne entre les trois équipes sur une semaine de travail
(hypothèse que les moyennes sont égales).
Nous allons donc ici appliquer la théorie vue en cours et vérifier si les résultats sont les
mêmes que dans Minitab et MS Excel. Nous utiliserons évidemment les mêmes données:
Nous pouvons d'abord faire une analyse qualitative avec le package ggplot2:
Nous avons alors d'abord en utilisant la commande native aov( ) native à R pour avoir le
modèle sans interactions:
Pour le modèle avec interactions (évidemment impossible dans le cas présent pour un logiciel
de communiquer les statistiques puisqu'il n'y a pas de répétitions):
et respectivement:
ou pour revenir avec aux machines avec un visuel plus proche de celui de Minitab mais aussi
meilleur que ce dernier:
Imaginons une entreprise faisant les trois huit. Nous avons trois équipes qui travaillent sur une
même machine. Nous souhaitons vérifier avec un seuil de confiance de 95% s'il y a une
différence de productivité moyenne entre les trois équipes sur une semaine de travail
(hypothèse que les moyennes sont égales).
Nous allons donc ici appliquer la théorie vue en cours et vérifier si les résultats sont les
mêmes que dans Minitab et MS Excel. Nous utiliserons évidemment les mêmes données:
Nous pouvons d'abord faire une analyse qualitative avec le package ggplot2:
Nous avons alors d'abord en utilisant la commande native aov( ) native à R pour avoir le
modèle sans interactions:
Remarque: Nous verrons juste après un test équivalent considéré comme plus robuste qu'est le
test de (l'étendue) Tukey HSD.
Nous partons des mêmes données que dans le cours théorique et que l'ANOVA à 1 facteur
fixe faite plus haut:
Le résultat peut surprendre dans le cas présent mais c'est simplement que la différence n'est
jamais significative et donc que les p-value sont in extenso très proches de l'unité.
Ici, nous allons vérifier si les concepts démontrés dans le cours de statistique théorique
concernant le test de Tukey sont appliqués de façon identique dans R par rapport à Minitab.
Nous n'avons pas en classe dans le cours théorique simulé les distributions de l'étendue
studentisée plus par flemme que par manque de temps (les simulations de Monte-Carlo se
basant toujours sur le même principe...), donc il n'y aura pas de comparaison par rapport avec
MS Excel. Nous ferons ici juste une vérification des résultats renvoyés par Minitab en
utilisant les tables des étendues studentisées.
Nous partons des mêmes données que dans le cours Minitab et que l'ANOVA à 1 facteur fixe
faite plus haut:
Nous retrouvons bien les mêmes valeurs que dans le cours Minitab avec un visuel plus sobre
et donc moins confus et donc aussi avec les mêmes conclusions.
Qu'il n'y a priori pas d'option (contrairement à Minitab) pour faire un test par paire
individuelles avec le taux d'erreur individuel de Fisher.
L'objectif ici va être le cas plus intéressant d'application du test de Levene aux données d'une
ANOVA empilée et de vérifier encore une fois si les résultats sont conformes aux
simulations de Monte-Carlo effectuées dans le cours théoriques et à l'application numérique
effectuée dans MS Excel et Minitab.
Nous allons donc ici réutiliser le test de Levene testLevene( ) du package car vu plus haut
mais donc aussi le test Bartlett qui sera nouveau dans le cas présent.
Nous retrouvons donc les mêmes valeurs qu'avec Minitab (et donc avec les mêmes
conclusions) ou que dans le cours théorique.
Même si nous n'avons pas fait la démonstration du test de Bartlett en cours puisque les détails
de celle-ci est introuvable dans les livres et même dans l'article d'origine de Bartlett lui-même
(donc si quelqu'un à la démonstration détaillée qu'il se manifeste!) regardons quand même si
nous retrouvons la même valeur qu'avec Minitab en utilisant la fonction bartlett.test( ):
Donc bien que ce soit beaucoup plus précis c'est le même résultat que celui donné par
Minitab.
Non sans mal nous avons démontré dans le cours théorique le concept et les développements
relatifs à l'ANOVA hiérarchique complète (full hierarchical ANOVA). Nous allons vérifier ici
que nous obtenons les mêmes résultats que Minitat et surtout du NIST lui-même.
http://www.itl.nist.gov/div898/handbook/ppc/section2/ppc233.htm
Machines
1 2 3 4 5
0.125 0.118 0.123 0.126 0.118
Operator 0.127 0.122 0.125 0.128 0.129
Day 0.125 0.120 0.125 0.126 0.127
0.126 0.124 0.124 0.127 0.120
0.128 0.119 0.126 0.129 0.121
0.124 0.116 0.122 0.126 0.125
0.128 0.125 0.121 0.129 0.123
Operator
0.127 0.119 0.124 0.125 0.114
Night
0.126 0.125 0.126 0.130 0.124
0.129 0.120 0.125 0.124 0.117
Nous pouvons voir que R ne renvoie pas la ligne des totaux ni le deuxième test de Fisher mais
sinon le reste de sortie correspond bien à la documentation du NIST:
Pour combler cette lacune nous allons alors procéder en deux étapes. La première ci-dessus a
nous a permis d'obtenir le premier test de Fisher, celle ci-dessous nous permet d'obtenir la
deuxième:
Sinon la conclusion est la même qu'avec Minitab car ce dernier nous affichait les p-value pour
les deux tests de Fisher: Les Machines sont le seul facteur significatif donc les seules
améliorations doivent être procédées sur l'outil et non sur l'humain.
R ne propose pas explicitement au même titre que Minitab de type d'ANOVA qui dans le
cadre des plans d'expériences pourrait être utilisé pour chercher la combinaison optimale…
donc il va être difficile de l'optimiser… La seule chose que nous pouvons faire a priori c'est
une analyse de la variance d'où la présence de ce sujet dans le chapitre des ANOVA.
Nous allons comme à l'habitude reproduire l'exemple du cours théorique et vérifier que nous
obtenons bien les mêmes résultats que ceux calculés à la main et obtenus avec Minitab. Pour
cela, nous partons des données suivantes reproduisant le tableau carré latin du cours
théorique:
Considérons par exemple le cas d'une analyse de rapidité de 4 véhicules différents {1,2,3,4}
par 4 pilotes différents {I,II,III,IV} sur 4 routes différentes {A, B, D, C} et les valeurs
mesurées sont les consommations de carburant pour parcourir une distance fixe:
Véhicules
1 2 3 4
I A B D C
19 24 23 26
II D C A B
Pilotes
23 24 19 30
III B D C A
15 14 15 16
IV C A B D
19 18 19 16
Dans R nous importons un fichier CSV structuré comme le montre la capture d'écran ci-
dessous:
Remarquons que nous pouvons reconstruire le carré latin (bon ici il est transposé mais cela ne
change rien sur le fond):
Nous pouvons déjà constater ci-dessous que les pilotes on une influence certaine alors que
pour les autres rien n'est moins sûr ce que nous verifierons avec les tests de Fisher:
Nous retrouvons bien les valeurs calculés à la main! Avec les mêmes conclusions: Seuls les
pilotes influencent significativement le temps écoulé!
Notre but ici est de vérifier si nous réobtenons bien les résultats des calculs effectués à la
main dans le cours théorique et ce avec quel niveau de détails.
D'abord, nous reprenons les données du cours théorique qui sous forme esthétique est donnée
pour rappel par:
Méthode A Méthode B
Sujet Xa Ya Sujets Xb Yb
a1 5 20 b1 7 19
a2 10 23 b2 12 26
a3 12 30 b3 27 33
a4 9 25 b4 24 35
a5 23 34 b5 18 30
a6 21 40 b6 22 31
a7 14 27 b7 26 34
a8 18 38 b8 21 28
a9 6 24 b9 14 23
a10 13 31 b10 9 22
Moyenne 13.1 29.2 18.0 28.1
C'est un peu basique comme résumé puisqu'il n'y aucun test post-hoc mais bon sinon nous
retrouvons bien les résultats calculés dans le cours théorique à main avec la même conclusion!
Nous allons à nouveau contrôler que nous obtenons bien certain éléments calculés à la main
dans le cours théorique.
Ensuite, nous calculons les moyennes de chacun de variables dépendantes pour chacun des
groupes:
Nous pouvons déjà constater qu'entre les varibles dépendantes ET entre les variables
indépendantes il y dans les deux cas de variations qui semble être très fortes.
Donc conformément à ce que nous pouvions nous attendre, nous rejetons l'hypothèse nulle.
Nous allons ici vérifier que nous retrouvons les résultats de Minitab et obtenus avec MS Excel
avec les relations obtenues lors des démonstrations mathématiques de la méthode dans le
cours théorique. Nous reprenons donc les données partielles des Iris de Fisher (histoire de
faire original...):
Donc à part la première ligne que nous n'avons pas calculée dans le cours théorique (et que
Minitab ne donne pas), le reste est parfaitement conforme aux calculs faits à la main et dans
MS Excel ainsi que dans Minitab.
7
Pour rappel il s'agit d'un "scree plot" (en anglais "scree" se sont les éboulis qui descendent des montagnes).
ou plus classiquement:
Ce qui est parfaitement conforme à ce que nous avons obtenu dans les autres cours.
Nous pouvons aussi obtenir la matrice des vecteurs propres (aussi parfaitement conforme):
Donc ici, contrairement à Minitab, nous retrouvons les valeurs telles que calculées dans le
cours théorique à la main pour les points par contre pour les vecteurs propres cela est différent
de ce que nous avons calculé dans le cours théorique ET différent de ce que Minitab ou
MATLAB retournent...
Le but va être ici de vérifier comme à l'habitude que les résultats calculés à la main (et donc
les conclusions y relatives) sont conformes aux démonstrations mathématiques faites dans le
cours théorique avec l'exemple particulier que nous avions choisi ainsi que de vérifier s'ils
sont conformes à ce que nous avions dans MS Excel et Minitab.
Nous utiliserons la commande fa( ) du package psych qui nécessite lui-même le chargement
du package GPArotation pour contrôler le type de rotation.
D'abord nous montrons les données que nous utilisons et la matrice de covariance y relative:
Nous retrouvons donc presque les mêmes résultats que dans le cours théorique et
Minitab/MS Excel. C'est principalement la communalité qui diffère et scores des facteurs à
quelques pourcents.
Par contre nous n'avons pas démontré ni parlé dans le cours théorique de toutes les
informations données en-dessous des deux tableaux ni de ce paramètre u2 que renvoie R et
qui semble (mais à vérifier...) être simplement 1-h2. De plus les noms choisis par R pour les
colonnes ne sont pas des plus subtils... Il en sera de même pour les deux autres méthodes que
nous allons voir ci-dessous et qui ont toutes deux décortiquées dans le cours théorique.
Donc encore une fois à quelques % près non significatif, nous pouvons considérer les résultats
identiques à ceux de Minitab et MS Excel ainsi qu'à ceux obtenus dans le cours théorique.
Nous avons survolé rapidement dans le cours théorique les notions de M-estimateurs
(estimateurs robustes de la moyenne). Comme il s'agit d'ingénierie statistique nous n'avons
pas insisté trop longtemps sur les mathématiques sous-jacentes mais surtout sur les idées de
leur construction et les contraintes.
L'idée est de vérifier que ces quelques estimateurs sont disponibles dans R (comme ils le sont
dans SPSS ou SAS) et de vérifier que nous trouvons les bons résultats:
Nous allons ici comparer le résultat donné par le package quantreg par rapport aux trois
situations calculées avec Microsoft Excel dans le cours théorique:
Rappelons que cette technique est très utile en finance et dans le domaine biomédical.
Cependant on lui préfère dans les cas non linéaires l'interpolation par spline. Voici par
exemple un cas typique de l'interpolation quantile par spline que nous avons presque tous
connus en étant petit quand nous allions chez le docteur pour mesurer notre taille en fonction
de notre âge:
Maintenant, nous faisons une régression linéaire quantile de la médiane la commande rq( ) du
package quantreg et regardons ce que nous obtenons par rapport à ce que nous avions obtenu
dans MS Excel pour la médiane:
Donc c'est pas mal du tout par rapport à MS Excel! Il en est de même pour d'autres quantiles!
Maintenant, nous faisons plot de la régression linéaire OLS paramétrique avec des régressions
linéaires quantiles pour différents quantiles:
Comme nous pouvons le constater et l'avions fait observer dans le cours de méthodes
numériques, des régressions quantiles peuvent donc bien se croiser!
Là encore nous allons vérifier la conformité des résultats avec toujours les mêmes logiciels et
les calculs effectués à la main suite à la démonstration mathématique faite dans le cours
théorique.
Nous obtenons donc la même chose que les calculs effectués avec Minitab, MS Excel et à la
main.
Nous voulons à nouveau vérifier les calculs effectués à la main dans le cours théorique (avec
l'aide de MS Excel) et comparer aussi les résultats avec Minitab.
Donc d'abord, R a comme Minitab pas la possibilité à ce jour d'imposer l'espérance de loi de
Poisson donc il va en calculer l'estimateur.
Donc d'abord nous importons le jeu de données histoire de pas l'écrire à la main:
Bon c'est bien joli mais ceci c'est conforme aux calculs faits à la main mais allons au point qui
nous intéresse:
Là encore nous allons vérifier la conformité des résultats avec toujours les mêmes logiciels et
les calculs effectués à la main suite à la démonstration mathématique faite dans le cours
théorique.
Les résultats sont donc conformes à la différence que la p-value est extrêmement précise dans
R...
Même si nous l'avons vu précédemment, nous allons toutefois refaire l'exemple vu dans le
cours théorique et que nous avions contrôlé à la main avec MS Excel et Minitab.
Nous prenons les données suivantes (les mêmes que dans le cours théorique) avec la
commande fisher.test( ):
Nous allons ici encore une fois recalculer l'exemple vu dans le cours théorique et que nous
avions contrôlé à la main avec MS Excel et Minitab.
Bon nous n'allons nous pas nous étendre sur cette technique empirique qui est sujet à débat et
que de plus nous n'avons pas validé avec les techniques de Monte-Carlo dans le cours
théorique et qui n'est pas disponible dans Minitab.
Rappelons que cette technique ne serait approximativement valable que pour des tableaux de
contingences 2 2 dont certaines observations sont inférieures aux 5 unités et l'ensemble des
observations inférieures à 20 unités (d'après les tests qui auraient été faits par plusieurs
personnes de la communauté des praticiens de la statistique).
Pour cet exemple nous allons reprendre la petite table utilisée pour le test de Fisher.
Évidemment c'est erroné de la faire puisque les différentes cases de ce tableau de contingence
sont dépendantes entre elles. Mais bon... c'est juste pour montrer que R va voir qu'il y a un
faible nombre d'observations et va donc in extenso appliquer la correction de Yates.
Nous avons donc deux experts qui ont notés des vins. Les notes sont rangées dans l'ordre
croissant des notes et selon le numéro d'étiquette d'identifiant du vin (bref le même exemple
que nous avons utilisé après avoir fait les démonstrations mathématiques de cet indicateur
dans le cours théorique).
Le résultat est juste mais c'est dommage qu'il n'y ait pas d'intervalle de confiance ni de p-
value. Cependant avec la commande Kendall( ) du package Kendall on comble ces lacunes
facilement:
Par contre pour la p-value nous sommes loin du compte (même si la conclusion reste la
même...).
Encore une fois, le but est d'appliquer les calculs et démonstrations faits à la main pendant le
cours théorique de méthodes numériques et de vérifier que nous retrouvons les mêmes valeurs
avec R (et comparer par la même occasion avec Minitab).
Comme à l'habitude, le but va être de vérifier ici que nous retombons sur les calculs faits à la
main dans le cours théorique suite à la démonstration mathématique de l'origine de l'alpha de
Cronbach.
Donc nous retrouvons bien les mêmes résultats que ceux calculés à la main et que Minitab.
Et nous n'obtenons pas du tout le même résultat que dans le cours théorique ou avec Minitab
mais c'est normal. Effectivement pour réobtenir la même chose, il nous faudra écrire:
et voilà!
Nous voulons appliquer le test des rangs signés de Wilcoxon pour 2 échantillons appariés
comme démontré dans le cours de statistique théorique bien que celui-ci ne soit pas disponible
explicitement dans Minitab et ce afin de comparer au résultat calculé avec MS Excel dans le
cadre de l'approximation par une loi Normale.
Nous utiliserons pour cela les mêmes données brutes que dans le cours théorique avec la
commande wilcox.test( ) et d'abord sans approximation:
Nous retrouvons donc exactement les mêmes valeurs que dans Minitab (messages
d'avertissement exceptés) mais qui comme nous le savons ne correspondant pas à ce que nous
avons calculé à la main. Donc faisons l'approximation:
Nous voyons que cela ne change rien. Certes la valeur V est la même que dans le cours
théorique mais toujours pas la p-value.
Nous utilisons la commande SIGN.test( ) du package BSDA où la valeur test par défaut du
paramètre md est nul si pas spécifié mais la seule chose qui nous intéresse ici est l'intervalle!:
Nous avons effectué deux séries de mesures avec deux méthodes différentes.
Nous souhaiterions savoir la différence sont significatives ou pas. Le but étant aussi de
contrôler si nous avons un résultat différent de celui calculé à la main en cours et dans
MS Excel.
Nous utilisons la commande SIGN.test( ) du package BSDA où la valeur test par défaut du
paramètre md est nul si pas spécifié:
Nous retrouvons toutes les valeurs vues dans le cours de statistique théorique et que dans le
cours Minitab. Nous rejetons donc au vu de la p-value (inférieure à 5%) que la différence n'est
pas significative.
Nous pouvons aussi construire l'équivalent à partir du test binomial exact mais en focalisant
uniquement sur la p-value du résultat renvoyé:
Le but du test de Mood n'est pas de contrôler si des données appariées peuvent être
considérées comme égales ou non en se basant sur leur médiane (test des signes), mais
simplement de vérifier sur la base d'un tableau de contingence du Khi-deux, si le nombre de
valeurs au-dessus ou en-dessous de la médiane de deux échantillons est significativement
différente et proviennent donc de deux populations différentes.
Attention!!! Dans R il existe un test appelé test.mood( ) mais cela n'a rien à voir avec le test
de Mood des médianes. Contrairement à Minitab, je n'ai pas trouvé de test de Mood des
médianes intégré à R ou à un package particulier à ce jour.
Nous allons voir un exemple et montrer que nous pouvons l'obtenir intégralement à partir du
test du signe et du Khi-deux.
Ce qui donne:
Et ensuite:
Nous voyons que bien que la conclusion soit la même, le test fait une approximation par la loi
Normale.
Encore une fois, nous souhaitons comparer les démonstrations faites dans le cours de
statistiques théorique avec l'application numérique approximée dans MS Excel et exact avec
Minitab pour comparer avec le résultat obtenu avec R.
Nous retrouvons bien la valeur V égale à 26 que Minitab ne nous donnait pas directement. La
p-value est sensiblement différente de celle de Minitab et de celle calculée dans le cours
théorique (et rappelons que celle calculée à la main dans cours théorique était aussi
sensiblement différente que celle donnée par Minitab). Contrairement à Minitab, nous avons
cependant un intervalle de confiance ce qui est bien pratique. Concernant l'estimation
ponctuelle de la médiane, nous retrouvons bien la même valeur.
Comme d'haaaabitudeeee (en chantant...) le but est de voir si nous retrouvons les résultats
numériques calculés avec MS Excel en utilisant les relations démontrées dans le cours de
statistique théorique.
Par contre nous n'allons par reprendre les mêmes valeurs qu'en dans le cours théorique car
l'exemple y était trop petit (pour des raisons pédagogiques). Nous allons encore une fois
utiliser la commande native wilcox.test( ) de R:
Nous voulons appliquer le test des rangs signés de Wilcoxon pour 2 échantillons appariés
comme démontré dans le cours de statistique théorique. Rappelons que celui-ci ne n'était pas
disponible explicitement dans Minitab mais que faisant plusieurs manipulations nous avons
tout de même pu l'obtenir et le comparer au résultat calculé avec MS Excel dans le cadre de
l'approximation par une loi Normale. Nous allons ici vérifier la correspondance avec R.
Pour cela, nous allons encore une fois utiliser la commande native wilcox.test( ) de R:
Donc à part l'intervalle de confiance que Minitab ne donnait pas, nous retrouvons les mêmes
valeurs que dans ce dernier. La seule différence réside toujours dans la p-value qui est
significativement différence que celle obtenue à la main dans le cours théorique.
Cas à 1 échantillon
Nous souhaitons ici simplement vérifier que nous obtenons les mêmes résultats que ceux
obtenus dans MS Excel avec la méthode Monte Carlo (et in extenso vérifier les résultats
obtenus lors de l'étude théorique du test de Kolmogorov-Smirnov) et Minitab.
En utilisant la commande ks.test( ) et avec les mêmes données brutes que dans le cours
théorique, nous avons pour le test de Kolmogorov-Smirnov dans le cas particulier d'un test de
Normalité:
Donc ici nous voyons un énorme avantage par rapport à Minitab déjà: nous pouvons choisir la
distribution selon nos désirs!!
Nous n'avons pas la même valeur de D que Minitab (qui était 0.212) par contre elle
correspond à celle calculée dans le cours théorique à la main (idem pour la p-value). En réalité
avec R il y a un petit piège pour retomber sur ce que fait Minitab, il faut écrire (dès lors cela
équivaut à la variante dite pour rappel de "test d'ajustement de Lilliefor"):
Nous pouvons faire une double vérification en utilisant le test de Lilliefor explicitement via le
package nortest et la commande lillie.test( ):
Bien que la valeur D soit la même, la p-value est très différente! Il faudrait investiguer pour
savoir d'où vient cette différence. Peut-être dans la manière de calculer l'écart-type??? A
suivre...
Cas à 2 échantillons
Nous avons mentionné dans le cours théorique l'extension triviale de pouvoir généraliser le
test de Kolmogorov-Smirnov à deux échantillons empiriques. Voyons si R propose cela (nous
n'avais pas fait de calcul à la main). Considérons le jeu de données suivant:
Dans le package MASS il y a une commande fitdistr( ) et qui va chercher avec les techniques
classiques d'optimisation non linéaire à maximiser la log-vraisemblance. Voyons cela avec les
quantités de notre fichier des ventes d'articles en utilisant un script qui exécutera les
optimisations de plusieurs distributions d'un seul coup:
Encore une fois nous allons vérifier que nous retrouvons bien le même résultat que ce que
nous avons obtenu dans le cours théorique lors des calculs à la main ou dans Minitab.
Nous retrouvons donc les mêmes valeurs que celles calculés dans le cours théorique ainsi que
dans Minitab et donc avec les mêmes conclusions!
Au niveau de la conclusion, la p-value étant beaucoup plus petite que les valeurs
traditionnelles critiques (10%, 5%, 1%) nous mettons donc en évidence le fait qu'il y a une
différence significative entre le groupe de contrôle et de test à travers les différentes strates.
Nous allons ici vérifier validité du test de Grubbs qui est implémenté dans R que nous avons
vu dans le cours théorique et dont nous avons appris à calculer les valeurs critiques à la main
pour n'importe quelle distribution. Nous allons aussi comparer le résultat à ce que renvoie
Minitab.
Nous retrouvons donc la même valeur de G que celles calculées à la main par Monte-Carlo
dans MS Excel et que dans Minitab (à partir de la version 17 pour ce dernier). Par contre pour
les p-value nous en sommes à des kilomètres (il y a 50% de différence).
Au fait c'est un piège car il faut ajouter un paramètre pour retrouves les valeurs que nous
connaissons:
Nous allons ici vérifier validité du test de Dixon qui est implémenté dans R que nous avons
vu dans le cours théorique et dont nous avons appris à calculer les valeurs critiques à la main
pour n'importe quelle distribution. Nous allons aussi comparer le résultat à ce que renvoie
Minitab.
Nous retrouvons donc exactement les mêmes valeurs que celles calculées à la main par
Monte-Carlo dans MS Excel et que dans Minitab (à partir de la version 17 pour ce dernier).
Le point qui est dommage cependant c'est que contrairement à Minitab, nous ne pouvons pas
choisir le niveau de confiance (du moins à ce jour).
Nous allons nous baser sur l'exemple pratique fait aussi en cours pour appliquer la pseudo-
démonstration (guère convaincante) que nous avons étudiée concernant ce test. Nous allons
comparer le résultat obtenu aussi par rapport à Minitab.
Nous retrouvons donc exactement les mêmes résultats avec les mêmes conclusions que dans
le cours théorique.
Au même titre que dans le cours théorique et Minitab nous allons prendre par hommage
l'excellent exemple original de l'article de Kruskal-Wallis:
Donc à part le faire que nous avons moins d'informations que dans le cours théorique et
qu'avec Minitab, les résultats donnés sont en parfait adéquation.
Nous avons étudié dans le cours théorique les distributions bivariées de la loi Normale et
mentionné les différentes variantes des distributions de Student.
Le but ici va être de générer des variables aléatoires Normales bivariées dans un premier
temps et de calculer ensuite le meilleur ajustement.
Dans un premier temps nous faisons comme dans le cours théorique en simulant les
réalisations 10'000 variables aléatoires uniformes qui sont donc l'inverse d'une distribution
Normale centrée réduite bivariée dont le coefficient de corrélation est -0.5 avec la fonction
normalCopula( ) du package copula.
Voici une première représentation pas forcément très utile (opinion très personnelle):
Et à la page suivante une représentation beaucoup plus intéressante et remarquez que nous
retrouvons ce que nous avons construit à la main dans le cours théorique avec le tableur
Microsoft Excel:
Régressions
Ce chapitre est dédié aux différentes techniques de régressions. Rappelons que nous avons vu
dans le cours théorique qu'il y en un relativement grand nombre. Nous allons passer ici en
revue que celle dont nous avons étudiés et démontrés les développements mathématiques
pendant le cours théoriques.
Pour ce nous allons d'abord utiliser les mêmes données8 que dans le cours MS Excel et
Minitab (nous allons aussi passer sous silence la problématique des valeurs manquantes):
Sans package dans un premier temps, nous avons en utilisant la fonction lm( ):
8
Rappel!!! Normalement pour pouvoir appliquer le modèle Gaussien de régression linéaire il faut pour chaque
point d'abscisse plusieurs mesures en ordonnées (au moins 10).
Pour extraire le coefficient de corrélation seul il suffit d'écrire (nous verrons un peu plus loin
comment extraire d'autres informations):
Pour obtenir la somme des carrés des erreurs il suffit d'utiliser la fonction deviance( ):
Nous pouvons également directement échantilloner les données d'origines pour faire une
régression sur un sous-ensemble:
Et avec la package ggplot2 nous traçons d'abord la meilleure droite de régression selon la
méthode des moindres carrés avec l'intervalle de confiance à 95% en utilisant la commande
stat_smooth( ) du package ggplot2:
Nous retrouvons bien la même chose que dans le cours MS Excel et Minitab.
Ce qui est intéressant aussi et très important c'est de vérifier l'homoscédasticité des résidus.
Nous pouvons faire cela visuellement avec la commane suivante:
Bon évidemment... dans le cas présent comme nous n'avons qu'une seule mesure par point
d'abscisse il est presque impossible de dire quoi que ce soit de qualitativement sérieux!
Avec la commande predict( ), nous pouvons nous faire des projections avec intervalle de
confiance pour une ou plusieurs valeurs:
Nous pouvons faire mieux en associant intervalle de confiance et de prédiction comme nous
l'avions fait dans le cours Microsoft Excel et Minitab:
Ce qui donne:
Et donc nous retrouvons bien la même chose que dans les cours MS Excel et Minitab.
Nous retrouvons la même chose que dans MS Excel et Minitab à la différence qu'il y a des
informations manquantes avec R qui sont intéressantes dans la pratique et que la mise en page
n'est pas des meilleures...
Nous pouvons cependant faire mieux pour récupérer les valeurs manquantes:
Nous pouvons aussi afficher un graphique un peu similaire à celui de Minitab mais plus
technique incluant le leverage (niveau de levier) et la distance de Cook:
Nous pouvons faire un test de l'homoscédasticité sur la régression des rédisus puisque si ces
derniers satisfont les hypothèses, la droite passant par les résidus devrait avoir une pente nulle
(ou non significative) ainsi qu'une ordonnée à l'origine nulle. Or, dans le cas présent:
Nous pouvons aussi obtenir directement l'effet de levier avec la commande hatvalues( ):
Nous avons étudié dans le cours théorique quelques M-estimateurs avec un regard très
critique sur le choix de ces derniers. Voyons comme il est possible avec R de faire une
régression robuste basée par défaut sur le M-estimateur de Tukey à l'aide de la fonction rlm( )
du package MASS en reprenant les mêmes données que l'exemple précédent:
à comparer avec le modèle non robuste où nous avions obtenu pour rappel…:
Dans la littérature (la bonne de qualité!), nous voyons souvent les régressions linéaires
univariées représentées avec la loi Normale à la verticale. Comme c'est une fonctionnalité
souvent demandée voici typiquement le code avec l'aimable autorisation de
Arthur CHARPENTIER pris de son billet de blog du mois d'octobre 2011:
http://freakonometrics.hypotheses.org/2280
Donc avec des données bidons cela donne d'abord au niveau des commandes.
#du classique:
regmdl=glm(Y~X,data=df,family=gaussian(link="identity"))
#on prépare x pour dessiner les lignes
x=seq(min(X),max(X),length=500)
lines(regLine,lwd=2)
upLine=trans3d(x,y1,rep(0,length(x)),mat)
lines(upLine,lty=2)
lowLine=trans3d(x,y2,rep(0,length(x)),mat)
lines(lowLine,lty=2)
C=trans3d(X,Y,rep(0,length(X)),mat)
points(C,pch=19,col="red")
Ce qui donne:
Et nous pouvons évidemment faire de même avec d'autres distribution (pensez au modèles de
régressions linéaires généralisées Binomiaux, ou Poissons).
Nous pouvons faire comme nous l'avons vu dans le cours théorie une régression uni ou
multivariée non linéaire par moindres carrées ordinaires.
Tiré de l'ouvrage Biostatistcal Analysis and Design using R de Murrey Logan voici un
excellent petit résumé des cas classiques des fonctions utilisées avec nls( ):
Là encore, nous allons vérifier que nous retrouvons les mêmes résultats que les calculs faits à
la main après la démonstration mathématique dans le cours théorique du test et que dans
Minitab et ce avec les mêmes données en utilisant la commande cor.test( ):
Nous retrouvons donc les mêmes valeurs à la différence que R donne des informations plus
pertinentes et plus précises que Minitab (du moins par défaut).
Nous allons partir ici des mêmes données que celles utilisées dans le cours théorique pour
encore une fois vérifier que nous retombons sur la même chose ou pas et le comparer au
résultat renvoyé par Tangra.
Maintenant faisons la même régression linéaire multiple que dans le cours théorique et
Minitab avec la commande lm( ) en utilisant les données ci-dessous:
La véracité de ce diagramme est facilement vérifiable avec MS Excel (mais dans un temps
beaucoup plus long...).
Nous pouvons évidemment aussi très facilement faire des prévisions à nouveau avec la
commande predict( ):
Encore une fois nous allons contrôler les calculs faits à la main dans le cours théorie et
obtenus aussi avec MS Excel et Minitab en utilisant toujours les mêmes données:
Nous retrouvons effectivement les valeurs calculées manuellement en cours et donc les
mêmes conclusions.
Ici nous allons encore une fois vérifier que nous retrouvons les mêmes résultats que dans le
cours Minitab et R et ainsi qu'avec les calculs faits à la main suite à la présentation de cette
technique empirique dans le cours théorique.
Nous voyons que nous avons moins d'options que dans Minitab ce qui est bien dommage car
les options proposées par ce dernier sont très pertinentes. Cependant dans R, nous avons le
AIC (Akaike Information Criterion) et que nous n'avons pas à ce jour dans Minitab ce qui est
bien dommage.
Nous voyons cependant que le meilleur modèle résultant est le même car R s'arrêt dès que le
modèle se dégrade (le but étant de minimiser l'AIC d'où le nom de la commande...).
p ( x)
DS ( p, q) p( x) log p( x) log p( x) p( x) log q( x)
q ( x)
et nous remarquons de suite que si p( x) q( x) la distance entropique est alors nulle. Akaike
s'est concentré sur le deuxième terme uniquement car il contient la distribution supposée:
p( x) log q( x)
Ce dernier est par construction toujours négatif puisque nous prenons le log d'une valeur
toujours inférieur à l'unité! Ensuite après quelques manipulations mathématiques discutables
(probabilités conditionnelles, approximations de Taylor, convergence en loi, etc.) nous
arrivons à:
AIC (k ) 2 L ˆ 2k
qui ne marche donc que pour des gros échantillons vu les hypothèses des développements
sous-jacents.
Dans tous les cas à sélection de modèles qui utilisent le même nombre de mesures k, il faut
donc privilégier le modèle qui a la plus petite valeur (donc si c'est négatif, plus c'est négatif
mieux c'est). Rappelons que le maximum de vraisemblance L décrit par construction d'autant
mieux un modèle que ce dernier est grand! Dans le cas de modèles qui n'utilisent pas le même
nombre de points de mesure k la sélection est à ma connaissance plus délicate à cause de la
présence du +2k…
Nous retrouvons donc bien le même résultat que dans Minitab et qu'à la main (ou que dans
MS Excel). Le seul point qui est un peu dommage c'est que R s'arrête à la première étape.
Cela aurait été bien d'en montrer plusieurs comme le fait Minitab.
Donc ici à l'aide des p-value on voit que c'est CoutC qui est le plus significatif, ensuite CoutB
et ainsi de suite. Nous voyons aussi en face l'AIC ce qui est toujours utile.
et ainsi de suite...
Là nous voyons que c'est CoutB qui est le plus significatif. C'est donc toujours conforme à
Minitab et aux calculs faits à la main. Nous pouvons donc ensuite, nous amuser à ajouter cette
variable au modèle:
et ainsi de suite...
Donc le meilleur modèle au sens du coefficient de corrélation ajusté est celui qui contient les
trois variables explicatives.
Bien que la conclusion soit la même que dans Minitab, le lecteur pour test que si nous ne
forçons pas l'ordonnée à l'origine comme étant nulle, la commande leaps( ) n'arrive pas à
s'exécuter (contrairement à Minitab où nous pouvons forcer ou non l'ordonnée à l'origine).
Comme quoi... comme on dit toujours: il n'y a pas de bons modèles. Il n'y a que des modèles
meilleurs que les autres!
Nous obtenons:
Nous poursuivons en faisant la même régression polynomiale que dans le cours théorique et
Minitab avec la commande lm( ):
Nous retrouvons là encore une fois les mêmes résultats qu'avec Minitab et le cours théorique.
La différence étant qu'avec Minitab nous avons automatiquement une série de graphiques
récapitulatifs du modèle.
Nous pouvons cependant compléter avec le package ggplot2 comme nous l'avons fait avant:
Nous pouvons revenir sur un cas plus compliqué comme nous l'avons fait avec Minitab pour
vérifier les possibilités de R. Nous allons donc utiliser les mêmes données que dans le cours
Minitab:
Le but ici va être de reprendre l'exemple fait dans le cours MS Excel utilisant la méthode de
Gauss-Newton et aussi la comparer avec ce que nous avions obtenu avec Minitab 16. Nous
partons donc des données suivantes:
Ensuite, nous utilisons la commande native pour y répliquer le modèle théorique supposé:
Nous remarquons donc que le nombre d'itérations est inférieur à celui de Minitab. Par contre
nous avons moins d'information concernant l'erreur du modèle que ce que donne Minitab ou
MS Excel de façon automatique. Avec un plus pour Minitab qui donne l'intervalle de
confiance au lieu de la valeur critique de la distribution de Student!
Sinon nous retrouvons exactement les mêmes valeurs de coefficients qu'avec Minitab et la
même erreur résiduelle.
Bon cet exercice n'est que pour les paillettes car avoir un modèle à seulement deux variables
explicatives est super limitatif (c'est toujours le problème des graphes car en dehors de la 2D
et de la 3D... c'est fini!).
Et maintenant nous jouons avec la commande scatter3d( ) et lm( ) des package car et rgl.
D'abord nous tentons de voir ce qu'un modèle purement linéaire avec deux variables
explicatives donne:
Soit graphiquement:
Et en mélangeant les deux (par contre à ce moment là on ne peut plus avoir les projections des
points sur le modèle):
Bref c'est joli mais très limite puisque à trois dimensions et de toute façon c'est aussi peu
précis (un bon tableau avec des chiffres en statistisques ça reste toujours la base!).
Comme nous l'avons vu dans le cours de statistique théorique, le problème d'une régression
linéaire simple (bivariée) ou multiple à variable expliquée de type binaire c'est que le modèle
théorique peut prendre toute valeur dans ce qui est évidemment une aberration puisque
dans le cas d'une variable expliquée de type 0/1, le modèle théorique ne doit jamais donner
une valeur inférieure à 0 et supérieure 1. Or c'est malheureusement ce que fera une régression
linéaire…
Pour cela nous avons démontré qu'il était possible d'utiliser une autre technique mathématique
qui nous assurait que la condition susmentionnée soit satisfaite. Nous allons reprendre le
même exemple mais au lieu de faire les calculs à la main, nous les faisons donc avec R (et les
comparer accessoirement avec R).
Nous utilisons la commande native glm( ) et nous remarquons que nous obtenons les mêmes
résultats que dans Minitab (avec la même différence que dans le cours théorique et pour les
mêmes raisons):
Nous pouvons construire un petit data frame avec les données du modèle:
Nous pouvons faire un plot des points mesurés, des points abscisses mesurées avec la théorie
ainsi que d'un lissage:
Nous pouvons également avoir la courbe ROC à l'aide de la commande roc( ) du package
pROC( ):
La courbe ROC est identique à Minitab et à celle calculée dans le cours théorique à la
différence que pour Minitab nous avons du laborieusement aller chercher une macro sur le
site web de l'éditeur du logiciel.
Avec la commande coords( ) du package pROC nous pouvons obtenir le meilleur cut-off:
Nous obtenons donc la même chose qu'avec Tanagra et qu'avec les calculs effectués à la main
dans le cours théorique.
Nous allons vérifier dans un premier temps si nous obtenons les mêmes résultats que ceux
calculés à la main et dans un deuxième temps comparer l'odds ratio et son intervalle de
confiance avec ce qui était donné par Minitab (ce dernier ne denant aucune information sur le
risque relatif) et le risque relatif et son intervalle de confiance à ce qui était donné par
MedCalc.
Nous importons donc d'abord les mêmes données que dans le cours théorique, nous les
mettons sous forme de table de contingence et en faisons un petit barplot( ) classique:
Rappel: Un odds ratio de 3.30 signifie pour rappel dans le cas présent que ceux traités au
placebo présentent un facteur de 3.30 d'exposition supplémentaire au risque cardiaque par
rapport à ceux traités à l'aspirine.
Nous avons donc avec ce package un odds ratio de 3.30 au lieu du 3.32 tel que calculé à la
main et avec Minitab et un intervalle de confiance de [1.48,7.89] au lieu de [1.55,7.14] tel que
calculé à la main et avec Minitab. La différence est donc raisonnable...
Pour le Inc risk ratio qui n'est en fait que el risque relatif, nous avons donc une valeur de 1.5
(au lieu du 2.21 tel que calculé à la main et avec MedCalc) et un intervalle de confiance de
[1.19,1.89] (au lieu de [1.266,3.869] tel que calculé à la main et avec MedCalc). Par contre
ici la différence est très suprenante! Si quelqu'un en connaît la raison...
Nous allons encore une fois vérifier ici si les calculs effectués à la main et dans MS Excel
suite aux démonstrations mathématiques d'une partie des éléments de la régression de Deming
effectuées dans le cours théorique.
Pour cela, nous utiliserons les 39 lignes suivantes (prises en partie de l'aide Minitab) qui
correspondent à la mesure de la pression sanguine avec deux instruments de 39 individus
(chaque individu ayant mesuré sa pression sanguine une fois avec l'ancien et une fois avec le
nouvel instrument):
Nous voyons alors que bien que la pente soit très proche de ce que nous avons obtenu avec
Minitab, au niveau de l'Intercept il y a près de 20% de différence!
Et là nous sommes beaucoup plus proches des résultats donnés par Minitab ce qui est heureux
et de plus, nous avons les intervalles de confiance avec donc les mêmes conclusions.
Nous pouvons aussi ploter les trois modèles de régressions (pour voir les différences quand
celles-ci sont éventuellement significatives):
Nous allons ici vérifier si nous retrouvons bien les différents résultats calculés à la main et
obtenus dans le cours théorique après avoir fait la démonstration détaillée de certains modèles
GLM particuliers (nous ne reviendrons donc pas sur le cas gaussien déjà traité plus haut).
Donc en ce qui concerne les coefficients nous obtenons les mêmes valeurs que dans le cours
théorique. Pour le reste... nous ferons confiance à R...
Nous nous retrouvons donc avec les mêmes résultats des coefficients et de l'ordonnée à
l'origine que dans le cours théorique et qu'avec SAS. Après pour les autres informations...
c'est une autre histoire!
On peut faire un plot du modèle. Voyons cela (si vous avez des propositions d'amélioration,
n'hésitez pas à me contacter!). D'abord montrons que nous pouvons obtenir les valeurs du
moyennes d'absence du modèle théorique pour chaque genre pour les notes moyennes:
Ensuite, nous pouvons faire cela pour toute l'étendue par genre et par note avec un graphe:
Ce qui donnera:
Donc il est rassurant de voir que plus les absences sont faibles, meilleure est la note...
Le but va être ici de vérifier si nous obtenons ou pas les résultats des calculs vu dans le cours
théorique lors de la lecture du l'ouvrage de M. Tenenhaus9 sur la régression PLS univariée
(PLS1), c'est-à-dire la régression sur des variables explicatives corrélées avec une unique
variable à expliquer.
Nous utiliserons les mêmes données que les calculs effectués à la main et que dans le cours
Minitab:
9
Michel Tenenhaus, Régression PLS, Édition Technip, ISBN 2-7108-0735-1, Pages 75-83
Nous allons demander à voir celles qui correspondent aux calculs faits à la main dans le cours
théorique et par Minitab:
Nous allons reprendre ici exactement le même exemple à hasard proportionnel non stratifié
dépendant du temps que celui fait dans le cours théorique après la démonstration
mathématique (pénible) des notions fondamentales de ce modèle....
Nous retrouvons donc le même résultat qu'à la main (et que pour MedCalc) pour le
coefficient. Pour le reste, cela sort à ce jour du contenu du cours théorique!
Interpolations
Les techniques d'interpolation sont contraitement aux techniques de régression des outils
mathématiques permettant d'ajuster au mois un courbe selon plusieurs critères empiriques
(comme pour les régressions) mais dont le but n'est pas ensuite de pouvoir faire des
projections (extrapolations) en avant et en arrière.
Les cas les plus courants d'interpolation sont les splines qui sont utilisées pour les courbes de
taux en finance. Voyons donc un exemple de cela.
http://www.federalreserve.gov/datadownload/Choose.aspx?rel=H.15
Nous avons alors avec beaucoup d'astuce d'abord en prenant la dernière date:
Et nous allons sur la dernière ligne droite en utilisant les commandes bs( ) ou ns( ) qui sont
différentes modèles de splines disponibles dans le package splines:
Il n'est pas possible (et cela n'a pas de sens mathématiquement parlant à ma connaissance) de
faire des projections (extrapolations) avec des interpolations par splines au-delà du domaine
de définition d'origine.
Bon ici nous n'allons pas trop insister car c'est de l'ingénierie numérique et comme nous
l'avons vue dans le cours théorique de Méthodes Numériques où nous avons implémenté le
code VBA du NIST dans MS Excel il existe de nombreux algorithmes LOESS (LOcally
Weighted regrESSion) différents qui donnent des résultats sensiblement différents.
Nous allons donc ici juste nous concentrer sur la mise en application de la méthode et le
résultat visuel sans le commenter:
Sondages
Fiabilité/Survie
Exercice 207.: Analyse de données censurées selon modèle non
paramétrique de Kaplan-Meier
R 3.0.2
Le but ici va être à nouveau de vérifier si nous obtenons bien les mêmes résultats que ceux
calculés à la main dans le cours théorique lorsque nous avons étudié la démonstration de
l'estimateur de survie de Kaplan-Meier.
Nous partons du tableau vu en cours mais adapté aux exigences de Minitab et R (il n'est
vraiment pas aisé de deviner que c'est sous cette forme que les choses doivent être
représentées):
Ensuite, en utilisant le package Survival, nous utilisons la commande survfit( ) avec Surv( )
de la manière suivante:
le ~1 signifiant qu'il n'y a pas de facteur explicatif particulier supplémentaires (comme l'âge,
le sexe ou autre).
Ce qui donne graphiquement (je n'ai pas trouvé comme désactiver les deux intervalles de
confiance sachant que nous n'en avons pas fait la démonstration mathématique dans le cours
théorique):
à comparer avec Minitab (on appréciera donc dans R l'indication des données censurées qui
sont a priori absentes dans Minitab):
Bon les graphiques c'est bien joli mais nous voulons des chiffres! Nous écrivons alors:
et nous retrouvons tous les chiffres calculés à la main dans le cours théorique (excepté pour
les deux dernières colonnes puisque nous n'avons pas démontré mathématiquement comme
obtenir l'intervalle de confiance).
Et en écrivant simplement:
Nous allons ici vérifier si nous retrouvons les résultats des calculs faits dans Microsoft Excel
suit à la démonstration faite dans le cours de Statistique théorique.
Nous partons donc des mêmes données que celles pour l'analyse de Kaplan-Meier puisque ce
sont celles que nous avons aussi utilisées pour le calcul dans Microsoft Excel après avoir fait
le calcul à la main de l'estimateur de maximum de vraisemblance de la loi de Weibull à deux
paramètres avec données censurées à droite:
Donc pour le paramètre de forme nous obtenons 1.35 (à comparer avec Microsoft Excel où
nous avions obtenu 1.18 et Minitab où nous avions obtenu 2.50) et pour le paramètre d'échelle
33.76 (à comparer avec Microsoft Excel où nous avions obtenu 27.16 et Minitab où nous
avions obtenu 26.02).
http://cran.r-project.org/bin/windows/base/old/
http://www.foxgo.net/10/post/2011/10/weibull-toolkit-for-r.html
Une fois l'installation de R 2.12.2 effectuée ainsi que celle du fichier Zip:
Plans d'expérience
Le but ici va être de voir dans quelle mesure nous pouvons reproduire dans R les plans
d'expérience vu dans le cours théorique et comparer les résultats au logiciel Minitab que nous
avons utilisé dans le cours de Génie Industriel.
Remarque: À ce jour (début 2014), R est très loin des possibilités des logiciels comme JMP
ou Minitab pour les plans d'expérience mais rien ne dit que dans quelques années...
Nous allons ici générer un plan d'expérience complet pour pouvoir prendre nos mesures avec
la commande expand.grid( ):
Ensuite, n'oubliez pas que vous pouvez éditer directement les données dans R avec la
commande edit(monplan). Ce sera tout pour le besoin le plus élémentaire imaginable...
Dans cet exercice, le lecteur n'aura pas besoin d'ouvrir de fichiers car le seul objectif ici est de
vérifier que R propose les plans d'expérience de base que nous avons étudiés dans le cours
théorique de Génie Industriel sur le même sujet avant de passer à des exemples concrets.
Concernant le fait d'aller plus loin avec les plans d'expérience, vous avez en gros deux
alternatives:
Pour les raisons déjà citées au début de ce e-book, je vais privilégier la deuxième méthode.
Plans factoriels
Plan factoriel complet général
Nous allons commencer par générer le même plan factoriel complet général que dans le cours
théorique de Génie Industriel en utilisant donc la commande fac.design( ) du package
DoE.base:
Nous voyons déjà dans la fenêtre de sessions que nous avons le même générateur et les
mêmes alias que ceux vus en cours.
Avec la commande aliasprint( ) on peut avoir les alias (qui sont donc les mêmes que Minitab
et que dans le cours théorique):
Nous ferons comme dans le cours théorique et avec Minitab, nous allons prendre 5 facteurs et
choisir le complet, le résolution V et le résolution III (comme le rappelle le schéma suivant de
Minitab):
qui lors du clic sur un des cas permet de récupere le plan dans l'espace de travail R.
Nous voyons bien les 16 essais comme dans le cours théorique et Minitab par contre les alias
ne son pas communiqués.
Maintenant voyons le pla factoriel fractionnaire ¼ de résolution III (plan de screening) Il est
très intéressant de voir comment les alias vont être choisis:
Le tableau et les alias sont conformes à Minitab mais la manière dont ils sont indiqués est très
vague (contrairement à Minitab) et peut porter à confusion.
Nous observons donc que la fenêtre de session n'affiche pas les alias (ce qui est dommage
puisque les plans de Plackett-Burman sont des plans à résolution III comme nous l'avons dit
en cours). Il n'est pas non plus possible de les obtenir avec la command aliasprint( ).
Le plan obtenu ci-dessus n'est pas le même qu'avec Minitab et ne correspond donc pas à
l'algorithme de Plackett-Burman que nous avons vu dans le cours théorique. Donc à creuser...
Le plan obtenu est lui plutôt surprenant. Nous n'avons que 4 colonnes pour les facteurs de
base. Certes il n'est pas utile d'en avoir plus mais cela aurait été sympathique de pouvoir
choisir...
Voici si jamais la liste des plans de Taguchi tels qu'indiqués dans l'aide de la commande (juste
pour info et en comparaison à Minitab):
Donc contrairement à Minitab ce package donne la table de Taguchi complète ce qui est plus
confortable intellectuellement. Profitons donc pour comparer par rapport à ce que nous avions
dans le cours théorique en désactivant la randomization:
Ce qui à part l'ordre, est parfaitement conforme au tableau obtenu dans le cours et où l'écriture
est beaucoup plus logique et intuitive à mon goût que dans Minitab:
La solution déterminée avec la vitesse et pression seule avait donc été de:
Package FrF2
Et nous voyons que bien qu'ayant défini un plan factoriel complet général, les coefficients:
sont donnés par rapport à un modèle aux variables codées (-1,1). Je n'ai pas trouvé la
possibilité comme sur Minitab de pouvoir choisir si le modèle final devait être donné en
variables codées ou non.
Donc on ne pourra pas comparer le plan factoriel général complet au plan factoriel général au
plan factoriel complet à 2 niveaux puisque ce deuxième nous est imposé!
Ensuite, nous passons d'abord à une analyse qualitative avec les fonctions wirePlot( ) et
contourPlot( ):
et nous voyons que le fait qu'on ne puisse pas bouger les légendes, changer les axes, tourner
les graphes est une désavantage majeur par rapport à Minitab:
Nous allons comme dans le cours Minitat considérer les plans CC comme des black box avec
des valeurs en entrées et en sortie dont nous connaissons les hypothèses d'utilisation.
Voyons d'abord ce que nous pouvons obtenir avec la fonction rsmChoose( ) du package
qualityTools:
Pour cela nous commençons par utiliser la fonction rsmDesign( ) avec les paramètres ad hoc
pour retomber sur le même plan que dans le cours Minitab:
Nous pouvons déjà mettre maintenant les noms des facteurs et leurs niveaux réels en utilisant
les commandes names( ), highs( ) et lows( ):
Ensuite vient la partie graphique où nous pouvons constater que le résultat sont très similaire
à Minitab (difficile de dire s'ils sont parfaitement identiques).
Plans de mélanges
Plan de mélange centré du simplexe
Le but va être de vérifier ici (comme à l'habitude) les calculs appliqués relatifs aux
démonstrations mathématiques faites dans le cours théorique relativement aux plans de
mélanges avec points aux sommets et comparer aussi accessoirement les résultats obtenus
avec Minitab.
La construction du plan de mélange étant faite, faisons en quelques analyses avec d'abord un
résumé:
Ensuite, nous appliquons le classique modèle linéaire généralisé gaussien avec interactions:
Alors nous retrouvons bien les valeurs des coefficients calculés à la main mais chose
curieuse... R (comme Minitab) n'arrive pas à calculer l'erreur ni les intervalles de confiance et
la p-value alors que nous avons démontré qu'un tel calcul était possible (et nous l'avons fait!).
Raison pour laquelle... (encore une fois!) il est important de rappeler qu'il faut toujours
plusieurs progiciels de statistiques pour faire des études sérieuses et un peu plus sûres.
Nous pouvons voir que celle-ci est plus intéressant qu'avec Minitab où nous n'avions qu'un
résumé objectivement peu exploitable (Minitab faisant la somme des termes linéaires et la
somme des termes quadratiques dans le tableau de l'ANOVA).
Finissons par un plot (le rendu est meilleur que Minitab mais les finitions prennent à ce jour
plus de temps que dans Minitab) en utilisant respectivement les commandes contourPlot3( )
et wirePlot3( ) du package qualityTools (ce qui nous permet d'éviter d'expliciter les termes
d'interactions car il les prendra en compte automatique)
Nous avons donc beaucoup moins d'options que dans Minitab mais pour du gratuit c'est quand
même extraordinaire.
La majorité de ce chapitre (si ce n'est la totalité en réalité) se satellise autour du package qcc.
Nous avons étudié dans le cours de génie industriel la notion de rendement global combiné et
de DPU/DPMO dans le cadre de processus auxquels nous appliquons l'approche Six Sigma.
Voyons ce que le package Six Sigma nous propose avec la commande ss.ca.yield( ) dans le
cas d'un procédé de fabrication où sur un total de 1915 unités produites nous avons à la
respectivement aux trois premières étapes du processus 3, 5 et 12 éléments défectueux et nous
avons à chaque étape respectivement réinjecté 1, 2 et 4 éléments retouchés (corrigé):
Contrôlons la deuxième valeur qui est le taux d'efficacité avec retouches (FTY: First Time
Yield):
Voyons comment le package SixSigma calcule la troisième valeur qui est le rendement global
combiné (RTY: Rolled Troughput Yield):
La définition utilisée n'est donc pas commune... puisqu'à chaque étape il reprend la quantité
initiale. Dommage... (du moins à ma connaissance).
Mais là encore... la définition prise par le package SixSigma n'est pas forcément des plus
judicieuses puisque normalement le DPU est le ratio du nombre total d'éléments défectueux
sur le nombre total d'unité de production observées...
Et pour le DPMO:
Raison pour laquelle il est toujours important de rappeler dans les rapports la définition
utilisée dans votre entreprise!
Voir par exemple (outre le support du cours théorique), le site suivant pour un définition assez
courante:
http://www.isixsigma.com/dictionary/dpo/
Nous allons vérifier ici que nous retrouvons donc les mêmes valeurs que celles calculées à la
main dans le cours théorique de Génie Industriel.
Rappelons le le contexte:
Dans un projet, une tâche de 4 jours a une marge totale de 2 jours. Si nous sortons de cette
marge, le projet va prendre du retard pour un coût de 1200.-/jour.
Quel est le coût réel associé à 1, 2 et 3 jours de retard selon le modèle de Taguchi?
Maintenant voyons un cas de calcul avec la situation vue dans le cours théorique, appelée "le
nominal est le meilleur.
Ensuite, nous utilisons la commande ss.lfa( ) du package SixSigma ce qui nous donne la
sortie suivante dans un premier temps:
Ce qui correspond bien aux calculs faits à la main dans le cours théorique. Nous avons aussi
automatiquement le graphique suivant donné par R qui correspond aussi à ce qui a été vu dans
le cours théorique:
Encore une fois, nous allons contrôler que nous obtenons les mêmes données que celles
calculées manuellement dans le cours théorique, ainsi qu'avec MS Excel et Minitab.
La fonction que nous allons utiliser ici n'est plus disponible depuis la version 2.2 du package
qcc (on se demande bien sous la pression de qui...). Il faut donc décompresser le fichier:
Ensuite, nous utilisons les commandes qcc.groups( ) pour préparer les données:
Ensuite, nous sommes obligés de choisir au moins une carte de contrôle pour l'analyse six
pack. Nous prendrons alors une carte X-barre en utilisant la commande qcc( ):
Donc à part l'esthétique qui est relativement médiocre mais améliorable en changeant le code
source, nous retrouvons exactement les mêmes résultats que ceux calculés dans le cours
théorique et dans Minitab excepté pour le LCL et l'UCL où il y a une différence de quelques
pourcents.
Et c'est maintenant que vient la commande qui a disparue dans les nouvelles versions du
package. Nous allons donc utiliser la commande process.capability( ):
avec le graphique:
Là aussi les résultats ne sont pas tout à fait les mêmes que ceux calculés dans le cours
théorique et que ceux calculés à la main mais il y a toutefois un avantage très intéressant avec
R: nous avons les intervalles de confiance pour tous les indicateurs de capabilité :-)
Avec le graphique associé (la encore on est loin derrière la qualité esthétique de Minitab où de
la pertinence des informations que donne ce dernier sur le graphe SixPack:
Maintenant, voyons ce que nous obtenons avec la fonction pcr du package qualityTools:
et avec le graphique est le indicateurs (dont nous voyons que ces derniers ne respectent pas la
notation d'usage... d'où l'importance de suivre les normes pour ne pas se mélanger les
pinceaux):
À comparer avec le graphique obtenu avec Minitab les données calculées à la main dans le
cours théorique!
Lors de notre étude des statistiques de rangs nous avons démontré qu'une constante émergait
qui dépendant du nombre de rangs et qui nous était indispensable pour construite un certain
nombre de cartes de contrôle.
Plutôt que de faire référence à des tables, le package SixSigma contient une commande
ss.cc.getd2( ) qui nous permet de récuperer la valeur de la constante de Hartley pour un
certain nombre de rangs (taille d'échantillon). Nous allons vérifier ici dans un premier temps
que nous obtenons bien les valeurs obtenues dans le cours théorique pour la première
constante de Harley d 2 :
De même pour la constante c4 dont nous avons démontré l'origine dans le cours théorique de
Génie Industriel:
Dans le cadre de contrôle de lots nous avons obtenus les mesures suivantes:
Nous obtenons donc les mêmes valeurs et la même carte que dans le cours de statistique
théorique et MS Excel ainsi que Minitab.
Cependant avec un plus: R détecte les anomalies de séquences sous la dénomination Number
violating runs. Sympa!!!
Dans le cadre de contrôle de lots nous avons obtenus les mesures suivantes:
Nous obtenons donc les mêmes valeurs et la même carte que dans le cours de statistique
théorique et MS Excel ainsi que Minitab.
Cependant avec un plus: R détecte les anomalies de séquences sous la dénomination Number
violating runs.
Dans le cadre de contrôle de lots nous avons obtenu les mesures suivantes:
Nous obtenons donc les mêmes valeurs et la même carte que dans le cours de statistique
théorique et MS Excel ainsi que Minitab.
Cependant avec un plus: R détecte les anomalies de séquences sous la dénomination Number
violating runs.
Dans le cadre de contrôle de lots nous avons obtenu les mesures suivantes:
Nous obtenons donc les mêmes valeurs et la même carte que dans le cours de statistique
théorique et MS Excel ainsi que Minitab.
Cependant avec un plus: R détecte les anomalies de séquences sous la dénomination Number
violating runs.
Dans le cadre de contrôle d'une production sous contrôle statistique nous avons obtenu les
mesures suivantes (chaque ligne est pour rappel une journée de mesures):
Nous obtenons donc les mêmes valeurs et la même carte que dans le cours de statistique
théorique et MS Excel ainsi que Minitab.
Cependant avec un plus: R détecte les anomalies de séquences sous la dénomination Number
violating runs.
Courbe d'efficacité
Nous avons vu dans le cours théorique comment calculer l'erreur de type II d'une carte de
contrôle de type X-barre. Le but ici est de vérifier que nous retrouvons les mêmes résultats
sous les mêmes hypothèses.
Dans le cadre de contrôle d'une production sous contrôle statistique nous avons obtenu les
mesures suivantes (chaque ligne est pour rappel une journée de mesures):
Nous obtenons donc les mêmes valeurs et la même carte que dans le cours de statistique
théorique et MS Excel ainsi que Minitab.
Cependant avec un plus: R détecte les anomalies de séquences sous la dénomination Number
violating runs.
Malheureusement il n'est pas possible à ma connaissance de combiner dans une seule fenêtre
la CC X barre-S et la CC S barre-S et ce même avec la commande par(mfrow=c(2,2))
contrairement à Minitab...
Il ne semble pas possible non plus, contrairement à Minitab, d'avoir une abscisse basée sur
des dates (timbres).
Dans le cadre de contrôle d'une production sous contrôle statistique nous avons obtenu les
mesures suivantes (chaque ligne est pour rappel une journée de mesures):
Nous obtenons donc les mêmes valeurs et la même carte que dans le cours de statistique
théorique et MS Excel ainsi que Minitab.
Cependant avec un plus: R détecte les anomalies de séquences sous la dénomination Number
violating runs.
Malheureusement il n'est pas possible à ma connaissance de combiner dans une seule fenêtre
la CC X barre-R (n'existe pas à ce jour) et la CC R barre-R et ce même avec la commande
par(mfrow=c(2,2)) contrairement à Minitab...
Il ne semble pas possible non plus, contrairement à Minitab, d'avoir une abscisse basée sur
des dates (timbres).
Dans le cadre de contrôle d'une production sous contrôle statistique nous avons obtenu les
mesures suivantes (chaque ligne est pour rappel une journée de mesures):
Nous retrouvons donc quasiment les mêmes valeurs que dans le cours théorique de
statistiques et que le cours MS Excel et Minitab. Effectivement, nous y avions fait une
approximation en ce qui concerne les limites de contrôle (indépendantes du temps t) et la
valeur du point de départ d'où la petite différence.
Dans le cadre de contrôle d'une production sous contrôle statistique nous avons obtenu les
mesures suivantes (chaque ligne est pour rappel une journée de mesures):
Nous voyons alors plusieurs curiosités. Bien que le forme de la carte CUSUM soit conforme,
impossible de définir a priori un V-masque et il semblerait qu'une deuxième ligne de mesures
est représentée dans la carte de contrôle (Ligne d'origine inconnue!). Les limites ne sont pas
non plus conformes à ce que nous avons calculé dans le cours théorique ou même obtenu avec
Minitab (le package prend le modèle des limites fixes – que Minitab propose aussi – mais cela
serait toutefois mieux d'avoir le V-masque puisque c'est l'intérêt principal de cette carte selon
moi).
La commande ne demande in extenso pas non plus les valeur de h et (ou in extenso de K).
A priori il vaudrait mieux mettre cette carte de côté à ce jour... Donc affaire à suivre!
Nous reprenons les mêmes données que dans le cours théorique et que dans le cours Minitab:
L'allure de la carte est conforme à ce que nous avons fait dans MS Excel et Minitab.
Nous avons:
La limite supérieure qui est de loin inférieure à celle de 278.2 donnée par Minitab
mais elle correspond parfaitement à la valeur calculée à la main dans le cours
théorique
La valeur centrale de 39.23 qui est bien supérieure à celle de 28.3 calculée par Minitab
mais elle correspond parfaitement à la valeur calculée à la main dans le cours
théorique
Nous allons vérifier ici si nous retrouvons les mêmes valeurs que dans le cours théorique et
que Minitab. Nous partons donc des mêmes données:
Nous retrouvons la même matrice de variances-covariances, la même limite LCL mais par
contre la limite UCL diffère considérablement et en plus nous ne savons pas pour quel niveau
de confiance il est calculé?!
Ce qui est le même graphique que celui obtenu à la main à part la limite supérieure qui est très
loin de celle calculée à la main ou de Minitab!
Pour ceux qui veulent l'ellipse de contrôle correspondant, ils peuvent se rendre à la page 531.
Une étude R&R mesure l'erreur de précision en faisant mesurer par exemple une pièce par
plusieurs fois et par des personnes différentes.
Étant donné que cette pièce ne change pas de taille, une variation dans les résultats doit
représenter la répétabilité de la mesure et sa reproductibilité par différentes personnes.
Cette étude est très simple et basée sur l'ANOVA à deux facteurs avec répétition dont on peut
déduire les données de l'étude R&R à partir du tableau de l'ANOVA (TAV). Pour plus
d'exemples, le lecteur pour se référer à la norme ISO/TR 12888:2011.
Considérons donc les données suivantes qui correspondent à une ANOVA à deux facteurs
avec répétition (il s'agit des mêmes valeurs numériques que l'exercice de l'ANOVA à deux
facteurs sans répétition vu plus haut mais détourné).
Nous reprenons donc les données du cours Minitab pour vérifier la concordance:
10
Voir aussi la fonction gageRRDesign( ) du package qualityTools
Nous retrouvons donc les mêmes valeurs que dans le cours Minitab. Le graphique nous donne
quant à lui (attention!! parfois il faut exécuter la commande ss.rr( ) plusieurs fois pour que le
graphique s'affiche sans erreurs):
La lisibilité du graphique n'est pas des meilleures sur un petit écran mais à part cela, nous
retrouvons bien les mêmes informations que celles renvoyées par Minitab.
Le but comme ici est comme à l'habitude de vérifier les calculs obtenus à la main dans le
cours théorique suite à la démonstration mathématique de l'origine du coefficient de
corrélation intra-classe ICC1.
Nous pardons donc des mêmes données que dans le cours théorique:
Nous utilisons le package psych pour faire l'analyse avec la fonction ICC( ) et le package
reshape2 avec la fonction acast( ) pour restructurer les données:
Nous retrouvons bien pour le ICC1 celui du calcul fait à la main suite à la démonstration
mathématique faite dans le cours théorique. Pour les autres ainsi que les intervalles de
confiance il faudra attendre la rédaction d'une démonstration mathématique élégante pour
qu'elle soit présentée en classe.
Bon revoilà le sujet qui fait mal et dont nous avons pas mal parlé dans le cours théorique de
MSP bien que la mathématique (empirique) sous-jacente soit triviale.
Ce qui donne:
Bon ce n'est pas surprenant on retrouve à un tout petit pourcentage près les mêmes valeurs
qu'avec Minitab mais rappelez-vous une chose. Dans l'état actuel, trois logiciels (Minitab, R
et QIMacros) utilisent trois définitions différentes (à ce jour) de la capabilité de gabarit. Donc
le jour où il y aura une norme internationale qui mettra tout le monde d'accord on en
rediscutera.
Bon ce sujet est comme nous le savons délicat étant donné qu'il existe plusieurs normes
nationales qui sont utilisées à l'instar des normes internationales ISO… Nous avons déjà pu le
vérifier avec Minitab où nous avions un écart d'environ 3% pour la taille de l'échantillon entre
les calculs théoriques et le logiciel et un écart de 7% pour le nombre d'éléments définissant si
le lot est à rejeter ou non.
Alors d'abord rappelons juste que dans Minitab nous avions pour les paramètres d'entrées ci-
dessous (pour un plan par attributs!):
Avec le package AcceptanceSampling il n'est pas possible de trouver le plan avec la même
approche. Ce package exige que nous utilisions les probabilités d'acceptation et les
pourcentages de défaillance. C'est-à-dire que nous supposons que nous avons à l'avanc les
données suivantes:
où:
PRP signifie Producer Risk Point, le premier argument est le niveau de qualité auquel
nous voulons évaluer le plan (AQL) et le second la probabilité minimale d'acceptation
du plan.
CRP signifie Consumer Risk Point, le premier argument est le risque de qualité auquel
nous voulons évaluer le plan (RQL) et le second la probabilité maximale d'acceptation
du plan.
Nous retrouvons alors des valeurs "relativement" proches du plan proposé par Minitab:
Nous allons voir comment utiliser les méthodes de Monte-Carlo (méthodes pseudo/stratifiées
ou de Fauré), du Bootstrapping, du Jacknife et du Latin Hypercube que dont nous avons
étudié les bases théoriques dans le cours.
Exercices: Monte-Carlo
Comme nous l'avons mentionné dans le cours théorique, un l'aspect central des méthodes de
Monte-Carlo est surtout relatif à la génération de variables aléatoires indépendantes ou non.
Nous allons ici utiliser la puissance du C++ derrière R et la simplicité de sa syntaxe pour
reproduire en quelques secondes un exercice pour lequel il nous a fallu un peu plus de 45
minutes à mettre en place dans le tableur MS Excel:
Avec un premier petit graph à l'aide de nos connaissances cumulée sur les graphiques:
et une deuxième petit graph encore plus important dans ce cas business:
Ce qui donne:
et ainsi de suite avec tous les graphs que l'on veut en quelques secondes. Personnellement je
trouve donc R plus intéressant et plus puissant que n'importe que tableur et même plus rapide
que @Risk.
Ce qui donne:
Méthode d'acceptation/rejet
Comme nous l'avons vu dans le cours théorique il s'agit d'une méthode dont l'algorithme
permettant de générer des variables aléatoires d'une fonction de densité difficilement
inversible est très simple à mettre en place.
L'idée très académique mais néanmoins plus réaliste qu'avec l'exemple pris dans la formation
MATLAB est de générer un vecteur de variables aléatoires d'une distribution beta de
paramètres 6 et 3:
Le but ici va être de mettre en pratique la technique jacknife mentionnée dans le cours
théorique de Méthodes Numériques à une petite régression univariée.
Nous allons de suite voir avec cet exemple pourquoi la taille de l'échantillon a une important
plus grande avec le Jacknife qu'avec le Boostrap que nous verrons juste après:
Nous allons ici mettre en oeuvre avec un facilité extrême l'exemple d'inférence par méthode
brute que nous avions mis en place avec MS Excel.
Nous partons du même jeu de données et allons aussi faire de l'inférence sur la médiane en
utilisant les commandes boot( ) et boot.ci( ) du package boot
Donc comme nous pouvons le voir pour la partie percentile que nous avions calculé avec
MS Excel nous ne sommes pas loin puisque nous avions obtenu (7.00, 29.5).... mais pas sans
peine!
Décisionnel
Bon pour l'instant pour être honnête je n'ai pas trouvé grand-chose au niveau décisionnel dans
R donc ce chapitre a un contenu qui est négligable. Cependant et comme à l'habitude cela
permet de vérifier des résultats calculés à la main dans le cours théorique.
Bon il s'agit ici d'analyse une matrice de préférence à un niveau qui pour rappele est basé sur
troix choix d'un facteur A:
[Choix A] 9 8 7 6 5 4 3 2 1 2 3 4 5 6 7 8 9 [Choix B]
[Choix A] 9 8 7 6 5 4 3 2 1 2 3 4 5 6 7 8 9 [Choix C]
[Choix B] 9 8 7 6 5 4 3 2 1 2 3 4 5 6 7 8 9 [Choix C]
D'abord vérifions qu'il y a une faible inconsistance a priori en analysant avec R les trois
valeurs propres:
Nous retrouvons donc avec R les valeurs calculées approximativement avec un logiciel
comme MS Excel excepté pour l'indice de Koczkodaj que nous n'avons pas étudié dans le
cours théorique.
Attention!!! Rappelons que dans la pratique, nous avons séparons normalement toujours les
données initiales en trois jeux: 1) Le jeu d'entraînement 2) Le jeu de test 3) Le tout. Ainsi,
L'ensemble des fonctions de Data Mining de R permet de mettre en entrée un jeu
d'entraînement et ensuite de faire une prédiction avec un jeu de test pour créer une matrice de
confusion et de conclure...
Il convient d'utiliser avec toute la prudence nécessaire et l'esprit très critique ces techniques
empiriques car parfois les résultats peuvent mener à des aberrations (et ce d'autant plus que le
jeu de données est petit).
Le Data Mining est la technique qui fait appel à la quantité chaque fois que l'on veut
congédier la qualité...
Signalons avec de commencer un petit package visuel pour faire quelques analyses de
Data Mining pour ceux qui aiment les interfaces graphiques... Il s'agit du package rattle:
qui donne:
11
Rappelons que le Big Data n'est pas une méthode d'analyse mais un champ d'étudre relativement au stockage,
à la gestion et l'exploitation des données.
La base rigoureuse du data minin est de séparer le jeu de données en 3 ensembles qui sont
dans l'ordre:
Une manière simple de procéder et d'utiliser la fonction sample( ) que nous avons déjà étudié:
Nous allons maintenant appliquer l'algorithme élémentaire One Rule rappelons que cet
algorithme fonction avec des variables (qualitatives) nominales et choisir parmi l'ensemble de
prédicteurs proposés seulement l'unique prédicteur qui en proportion explique le mieux la
variable d'intérêt.
Comme nous pouvons le constater (c'est tellement simple que même un enfant pourrait mettre
cet algorithme dans n'importe quel tableur) entre Déplacements et Moustaches, l'unique
Bref c'est sympa et ludique mais tellement trivial que nous n'irons pas plus loin.
Nous allons ici vérifier la technique de clustering ID-3 que nous avons étudié dans le cours
théorique de Méthodes Numériques et calculé à la main.
Nous allons donc travailler avec le fichier suivant et donc avec les mêmes données que dans
le cours théorique:
Nous n'avons donc pas a priori le niveau de finisse de paramétrage de Tanagra mais bon…
Dommage (au même titre que pour Tanagra) qu'il n'y ait pas de diagramme cependant... cela
aiderait à la compréhension.
Nous allons donc vérifier ici si nous obtenons le même dendrogramme que dans le cours
théorique (où nous nous étions follement amusé à calculer à la main....) et par la même
occasion nous pourrons vérifier si Minitab utiliser avec l'OBS une agglomération naïve ou
non. Nous utilisons donc les mêmes données que dans le cours théorique:
Nous pouvons d'abord avec la commande native dist( ) afficher la matrice symétrique des
distances euclidiennes:
Bon nous n'y voyons pas toute la matrice (pas assez de place sur la capture d'écran) mais les
valeurs sont bien les mêmes que celles calculées à la main. Vérifions maintenant que nous
avons bien le même dendrogramme avec la commande hclust( ):
Des méthodes d'analyses instables peuvent voir leur pouvoir prédictif amélioré par
perturbations et combinaisons des différentes résultats et en incluant des corrections
successives via des pondérations et différentes approches empiriques de corrections. L'idée du
bagging (ou "bootstrap aggrégatif") est donc similaire au boostrapping à la différence que des
paramètres ajoutés au modèle sont changés à chaque itération.
Effectuons un bagging du CAH précédent en utilisant une méthode proposée par Simon Raper
sur le lien suivant:
https://github.com/SimonRaper/buster
Ce qui donne:
Avec:
Nous reprenons les données utilisées dans le cours théorique (dont nous avons enlevé
l'appartenance aux groupes catégoriels connus a priori pour des raisons évidentes de
simplicité):
Que nous souhaitons associer en trois groupes distincts. Nous utilisons alors la commande
kmeans( ):
Nous retrouvons donc les mêmes résultats qu'avec Minitab (et a peu près les même que ceux
calculés à la main dans le cours théorique).
Et avec un autre visuel un peu plus techniques (je dois encore déterminer ce que les nombres
tout à droite signifient...) en utilisant les commandes daisy( ) et silhouette( ) des packages
respectifs cluster et HSAUR:
Rappelons que ce modèle fonctionne particulièrement bien avec un très grand nombre de
données et surtout si les données utilisées pour la distance sont numériques et continues afin
de pouvoir calculer la distance euclidienne comme c'est le cas ci-dessous. Dans le cas de
valeurs nominales ou ordinales, il faut avec la plus grande prudence les convertir en facteurs
numériques.
Et nous souhaiterions savoir comme selon la distance euclidienne, une tomate de douceur 6 et
de croquant 0 sera donc classifié... Nous obtenons alors:
Nous voyons que dans tous les cas la tomate est classifiée en tant que fruit (ouf!) que son plus
proche voisin (modèle 1-nn) fait partie forcément que d'une seule des classes (d'où le
1=100%). Avec les deux plus proches voisins (modèle 2-nn) nous avons encore une fois deux
plus proches voisins faisant partie de la même classe (donc 2/2=1=100%). Avec les trois plus
proches voisins (modèle 3-nn) nous avons 2 éléments de la même classe et 1 de classe
opposée (donc 2/3=0.66=66%). Avec les 4 plus proches voisions (modèle 4-nn) nous avons 2
éléments de la même classe et deux éléments de deux autres classes (donc 2/4=0.5=50%).
Nous allons ici encore et encore une fois utiliser les mêmes données que dans le cours
théorique:
Attention!!! Le package e1071, ne s'installera pas dans R3.0.2 ou d'autres versions si vous ne
spécifiez pas le serveur source comme étant l'Autriche (Austria) qui est le pays d'origine de
ceux qui ont développé ce package.
Nous pouvons même faire un contrôle sur les données d'entraînement elles-mêmes:
Donc comme nous pouvons le voir... un modèle sur une technique empirique reste une
modèle et de plus avec aussi peu de données...
Nous allons ici encore et encore une fois utiliser les mêmes données que dans le cours
théorique:
Rappel: Les techniques de classification (ségmentation) par arbres sont souvent utilisées en
marketing pour analyser les cibles répondant le mieux à certaines campagnes marketing. Les
praticiens du marketing appellent cela des "modèles uplift".
Pour cela, nous allons utiliser la commande rpart( ) du package rpart avec des paramètres
qui comme nous allons le voir sont loin d'être triviaux à deviner.
Comme nous pouvons le voir, le résultat n'est vraiment pas intuitif si nous comparons avec ce
que nous avons fait dans le cours théorique. Comme c'est choquant, nous allons donner un
poids plus grand à la surface qu'aux revenus:
et là c'est parfait!!! Par contre l'arbre s'arrête bien trop vite. Nous allons donc chercher à
l'étendre:
Pour information, minbucket=1 signifie que l'algorithme s'arrête s'il y a moins d'une feuille
dans la branche. Si vous utilisez le paramètre minsplit=10 (pas présent dans l'exemple ci-
dessus) cela signifie que l'algorithme ne coupera pas une branche s'il y a moins de 10
éléments dedans.
Mais nous avons vu dans le cours théorique que nous pouvons creuser l'arbre encore plus.
Alors qu'en est-il avec ce package? Essayons:
Que céééé beau! Nous obtenons donc la même chose que dans le cours théorique et avec un
niveau de détail supérieur à celui de Tanagra.
Remarquons que la commande rpart( ) utiliser bien par défaut l'indicateur de Gini pour vu
dans le cours théorique séparer les données. Effectivement, si nous les spécifions:
Maintenant faisons du "tree post-pruning". C'est-à-dire, prendre l'arbre qui minime l'erreur de
catégorisation. Pour cela, nous allons d'abord lister les arbres possibles à l'aide de la
commande printcp( ):
Donc il y aurait deux arbres qui minimisent la xerror qui sont respectivement celui à 4 et 7
sectionnements. Si nous voulions voir celui à 4 (plutôt que de chercher par tâtonnement) il
suffit alors de spécifier un CP à l'algorithme compris entre 0.010000 et 0.08333. Voyons cela
en utilisant la commande prune( ):
et ainsi de suite...
On peut ensuite faire des prédictions de clustering avec des données test en utilisant la
commande predict( ):
Donc le modèle est très bon. Nous pouvons calculer le taux d'erreur assez rapidement:
Donc comme nous l'avons mentionné dans le cours théorique, les Random Forests (forêts
aléatoires, appelées aussi "bootstrap Forests") ne consistent qu'en un bootstrapping de la
population initiale et ne garde à la fin que le modèle ayant le meilleur pouvoir prédictif!
Le package randomForest bien que fonctionnant très bien ne représente pas les données sous
forme structurée du meilleur arbre. Pour cela, nous allons utiliser le script suivant:
#**************************
#Build tree structure for Random forest
#Author: Andrea DAL POZZOLO
#**************************
getConds<-function(tree){
#store all conditions into a list
conds<-list()
#start by the terminal nodes and find previous conditions
id.leafs<-which(tree$status==-1)
j<-0
for(i in id.leafs){
j<-j+1
prevConds<-prevCond(tree,i)
conds[[j]]<-prevConds$cond
while(prevConds$id>1){
prevConds<-prevCond(tree,prevConds$id)
conds[[j]]<-paste(conds[[j]]," & ",prevConds$cond)
if(prevConds$id==1){
conds[[j]]<-paste(conds[[j]]," => ",tree$prediction[i])
break()
}
}
}
return(conds)
}
prevCond<-function(tree,i){
if(i %in% tree$right_daughter){
id<-which(tree$right_daughter==i)
cond<-paste(tree$split_var[id],">",tree$split_point[id])
}
if(i %in% tree$left_daughter){
id<-which(tree$left_daughter==i)
cond<-paste(tree$split_var[id],"<",tree$split_point[id])
}
return(list(cond=cond,id=id))
}
Nous pouvons aussi vérifier que les randomForest performent très bien à la classification:
Nous avons parlé dans le cours théorique de cette méthode empirique (à la base conçue pour
ne fonctionner que pour des variables catégorielles) de partitionnement/classification. Il n'y
avait rien à démontrer mathématiquement car elle n'utilistise que des éléments statistiques
triviaux.
Voyons-en un exemple pratique avec le jeu de données suivant où puisqu'il y des variables
non catégorielles nous allons devoir fractionner les valeurs continues et catégories en nous
basant sur notre retour d'expérience (pour rappel il s'agit là de la faiblesse majeure de cette
méthode):
Nous téléchargons d'abord le fichier *.zip de ce package depuis le site web suivant car à ce
jour il ne peut pas être installé depuis l'interface directement:
Le résultat n'est guère convaincant avec ce package. Au fait nous pouvons deviner qu'il utilise
le package rpart et donc il s'arrête au premier niveau pour les mêmes raisons que celles vues à
la page 942 (SPSS fait de même):
Mais ici comme nous ne pouvons a priori définir des poids nous sommes coincés pour cet
exemple scolaire.
Dès lors, prenons un jeu de données par défaut du package CHAID et jouons avec:
Ce qui donne:
Nous n'avons malheureusement pas le temps dans le cours théorique de vérifier que le résultat
soit juste… La comparaison avec SPSS montre que le résultt test très différent entre les deux
logiciels… mais cela est du au fait que nous avons pris un échantillon de taille 1'000. Si nous
prenons tout le jeu de données cela nous donne:
Ce qui est illisible (contrairement à SPSS) mais qui correspond par contre au niveau des
nœuds à ce que donne SPSS dont voici un morceau à la page suivante:
L'analyse d'affinité est très utilisée sur les sites Internets implémentant des modèles naifs de
proposition d'articles par rapport aux achats précédents de consommateurs.
Voyons en quoi consiste l'idée en partant du fichier suivant où chaque ligne représente un
panier commandé par un client donné (certains clients sont revenus plusieurs fois pour des
paniers différents) et où les colonnes représentées les articles achetés et les composantes les
quantités correspondantes:
Une fois ceci fait, nous préparons la matrice d'affinité et nous renommons les colonnes au
passage pour que leurs intitulés soient plus petits:
Pour le confort de lecture de cette matrice de distance nous allons trier les colonnes de plus
influentes au moins influentes et en refaire une matrice de distance (mais cette petite
transformation est totalement facultative et n'est pas une nécessité pour les graphiques que
nous verrons plus loin):
Nous avons calculé "à la main" dans Microsoft Excel la matrice des distances et les vecteurs
de position plongé dans un cas particulier de 2 à partir d'une matrice de dissimilarité et le
but ici va être de vérifier si quelques fonctions diponibles dans R à cet effet donnent les
mêmes valeurs ou permettent au moins d'arriver aux même conclusions (sachant qu'il de
nombreuses méthodes empiriques d'implémentation).
et graphiquement:
Et visuellement:
Et maintenant voyons que les les coordonnées peuvent vraiment varier suivant la méthode
mais que la conclusion qualitative reste bien la même. L'avantage des fonctions du package
vegan que nous allons voir maintenant est que nous retrouvons la fameuse variable de Stress
explicitement avec l'optimisation relativement à sa minimisation (nous nous passerons de la
construction du graphique qui est la même qu'avec la première méthode):
Nous allons voir ici si nous retombons avec R sur les calculs effectués dans le cours théorique
où nous avons utilisé la LDA avec le discriminant de Fisher.
Nous utiliserons donc les mêmes données et pour des raisons pédagogiques nous allons
d'abord le faire sans package:
Nous retrouvons donc exactement la même valeur et le même vecteur propre que dans les
calculs effectués à la main dans le cours théorique.
Ce qui donne:
A priori nous ne retrouvons pas le vecteur propre mais en réalité il n'en est rien. C'est juste
qu'avec la méthode vue dans le cours théorique le vecteur propre était normalisé alors qu'ici
ce n'est pas le cas. Nous pouvons vérifier que le vecteur LD1 est bien colinéaire au vecteur
propre obtenu avec la méthode vue dans le cours théorique:
Maintenant voyons qu'il est possible d'utiliser la LDA pour faire de la classification (via la
distance à la moyenne des groupes) et obtenons par la même occasion une matrice de
confusion:
Avec la commande partimat( ) du package klaR nous pouvons afficher les polygones de
spéparation des groupements:
Ce qui donne la frontière (droite) qui maximise les distances entre les points des groupes:
Le but ici va être de comparer le modèle de neurones construit dans le cours théorique avec le
tableur Microsoft Excel (avec pour avantage que l'on peut construire un réseau de neurones
multicouches avec n'importe quelle fonction de cut-off) par rapport à celui que nous allons
obtenir ci-dessous avec la commande neuralnet( ) du package neuralnet (il y a d'autres
package comme nnet, RSNNS ou encore AMORE) en utilisant le même jeu de données:
Nous ne nous attarderons pas trop longtemps sur cette fonctionnalité dans R car les neurones
à choix (quel que soit le package) sont très basiques et le manque de flexibilité de
construction fait que même un tableur performe mieux (c'est d'ailleurs ce que nous allons
pouvoir constater ici).
avec le graphique associé (n'oubliez pas qu'à chaque fois que vous relancerez le modèle,
les poids ne seront pas les mêmes!):
La somme des carrés des erreurs est de 3.9602 donc plus élevé que ce que nous avons obtenu
avec le tableur. Faisons donc un réseau à deux couches:
Aucune amélioration significative... et nous pourrions continuer ainsi longtemps sans obtenir
aucune amélioration.
Regardons cependant si ce modèle performe mieux sur un jeu de test n'étant pas inclus dans le
jeu d'entraînement. Nous allons utiliser le même jeu que dans le cours théorique:
Comme nous pouvons le voir, ce réseau de neurone a un problème puisqu'il considère chaque
entrée comme une valeur binaire. Ainsi, s'il existe un nombre en entrée il met celle-ci à la
valeur 1 et ensuite il multiplie par le poids et fait la somme des biais, raison pour laquelle on
se retrouve toujours avec la même valeur en sorite dans le cas présent. En réalité ce package
est fait pour de la classification et moins pour prédire des valeurs des continues...
Nous allons voir si nous pouvons faire mieux avec un autre package ((n'oubliez pas qu'à
chaque fois que vous relancerez le modèle, les poids ne seront pas les mêmes!):
Ensuite, comparons déjà la valeur prédite aux valeurs originales pour le jeu d'entraînement:
Ce qui pour cette exécution est moins bon que ce que nous avons obtenu avec le tableur mais
si vous exécutez plusieurs fois vous verrez que la valeur peut diminuer en-dessous de 0.8 que
nous avions obtenu avec Microsoft Excel. Il est aussi possible d'augmenter le nombre d'unités
dans la couche neuronale.
Donc nous faisons juste un peu moins bien... mais c'est le hasard car si vous relancez vous
verres que vous ferrez mieux que le réseau de neurone créé dans le cours théorique.
Nous savons que R met tout dans la mémoire RAM par défaut, or si cette dernière est limitée
par rapport aux données à analyser, outre acheter de la RAM (ce qui est compliqué dans les
multinationales car prend du temps) une possibilité et d'utiliser des packages qui gèrent les
données différemment (sur le disque).
Les fonctions detect_dm_csv( ) et laf_open( ) du package LaF (Large ascii File) ont pour but
premièrement de détecter le type de données des colonnes en espérant que le fichier soit pas
trop mal structuré (c'est souvent les deux premières lignes qui déterminent le résultat final) et
en second de créer la connexion au fichier avec les types de données détectés.
Nous avons alors dans R (il vaut mieux ne pas travailler sur l'ordinateur pendant le
traitement):
Nous voyons que nous avons importé environ 4.5 millions de lignes (ce n'est pas énorme) et
227 colonnes (ce qui est conséquent pour un sondage).
Nous voyons qu'il s'agit d'un objet de type ffdf content un dataframe. Malheureusement, le
fichier d'origine ne semble pas être bien structuré car la colonne AGEP pour "Age Population"
est censée être de type Integer et tout a été identifié comme des chaînes de caractères.
Après le principe reste le même que ce que nous avons vu jusqu'à maintenant cependant on
tombe souvent sur des abérrations avec ces gros jeux de données... :
Il est effectivement difficilement imaginable qu'il y ait 50'771 personnes dont l'âge est zéro et
ayant répondu...
Lorsque l'on traite des Téraoctets de données il devient difficile de travail avec un unique
serveur. Une solution consiste alors à utiliser un cluster et pourquoi ne pas utiliser celui
d'Apache: Hadoop. Wikitrends est un exemple d'utilisation de Hadoop analysant en temps réel
plus de 13 Téraoctets de données. La page officielle est ici:
http://www.wikitrends.eu/#/
On utilise la commande ts( ) du package zoo pour générer rapidement une série de dates sur
une base périodique donnée (frequency=12 donnera des mois, frequency=4 donnera des
trimestres et ainsi de suite):
Il est possible d'extraire des sous-parties en spécifiant la date de début et de fin de la période
de la série qui intéresse avec la commande window( ) du package zoo:
On peut s'amuser à obtenir le jour où le prix de l'action a été le plus élevé avec les techniques
de filtre habituelles de R:
Nous pouvons rechercher le maximum dans une certaine fenêtre de temps et calculer
basiquement la VaR historique:
Nous allons ici reprendre le même exemple que dans le cours théorique en utilisant la
commande decompose( ) du package zoo:
Nous n'avons pas étudié dans le cours théorique comment décomposer une série avec un
modèle multiplicatif mais voyons quand même comment faire (avec exactement les mêmes
données):
Faisons à nouveau la même chose mais avec un cas trivial pour voir comment les algorithmes
de décomposition performent en additif:
Avec la commande ts.plot( ) du package zoo, nous pouvons plotter que les composantes et les
mélanges qui nous intéressent:
Comme à l'habitude le but va être de vérifier que le logiciel utilise les résultats démontrés
dans le cours théorique avec un exemple implicitement lié à une série temporelle (les "0" étant
des "baisses" et les "1" étant des valeurs "haussières"). La question étant de savoir si les
séquences peuvent être considérées comme aléatoires ou non statistiquement parlant.
Nous devrions aussi retrouver les calculs faits à la main et à l'aide de MS Excel.
Bien que nous retrouvions les calculs faits dans le cours théorique, il est dommage que R ne
renvoie pas plus d'informations... Je préfère Minitab à ce niveau là et ce même si la
conclusion reste la même!
Nous allons ici vérifier que nous obtenons bien les mêmes résultats que ceux obtenus dans le
cours théorique et donc avec les mêmes données brutes et les comparer accessoirement avec
ce que nous avons obtenu dans MS Excel et Minitab.
Nous allons utiliser pour cela la commande filter( ) du package zoo et les mêmes données que
dans le cours théorique (données visibles ci-dessous dans la colonne data)
Donc tout est parfaitement conforme. Au niveau graphique cela donne donc:
Nous utilisons les mêmes données que dans le cours théorique nous laissons R trouver la
meilleure valeur du paramètre du modèle en utilisant la commande HoltWinters( ) native de
R:
Nous retrouvons donc la même valeur du paramètre alpha que dans le cours théorique (où
nous avions fait les calculs avec MS Excel) et dans le cours Minitab.
Pour la valeur de 256.2158 à quelques pourcents près non significatifs c'est à peu près
conforme au cous théorique et à Minitab. Mais attention à ces petites différences lorsque l'on
gère des projets dont les ordres de grandeur est le milliard de $...
Par contre il est dommage que R ne nous indique par les valeurs classiques d'erreurs que sont
les MAD, MSD et MAPE. Mais nous pouvons retrouver le MSD avec le calcul simple
suivant:
Donc à un pourcentage non significatif près, nous retrouvons la même que dans le cours
théorique et qu'avec Minitab.
Nous pouvons utiliser la commande native predict( ) pour faire des projections. Ici nous
allons faire un exemple aberrant allant au-delà d'une période:
Nous utilisons les mêmes données que dans le cours théorique nous laissons R trouver la
meilleure valeur des paramètres du modèle en utilisant la commande HoltWinters( ) native
de R:
Ainsi les deux coefficients sont très différents (+ de 40%) de ceux calculés à la main dans le
cours théorique. Par rapport à Minitab, la différence est de l'ordre de 5%.
Évidemment le meilleur choix est celui qui minimise le SSE. Ainsi dans le cours théorique,
nous avions obtenu un SSE d'un peu plus de 10'000... et avec Minitab de l'ordre de 6430...
donc R fait mieux si l'on peut dire...
Je n'ai par contre pas trouvé comment faire un lissage exponentiel double selon la méthode de
Brown.
Nous pouvons utiliser la commande native predict( ) pour faire des projections. Ici nous
allons faire un exemple aberrant allant au-delà d'une période:
Nous utilisons les mêmes données que dans le cours théorique et nous laissons R trouver la
meilleure valeur des paramètres du modèle toujours en utilisant la commande HoltWinters( )
native de R:
Nous ne pouvons pas ici comparer avec Minitab puisque ce dernier n'est pas capable de
trouver les 3 paramètres optimaux. Par contre par rapport faits à la main dans le cours
théorique il n'y a aucune comparaison possible... relativement différent pour les paramètres de
lissage où nous avions obtenu dans le cours théorique:
0.329, 1.06, 1
et:
Mais comme nous l'avons mentionné dans le cours théorique, ce modèle est aussi sensible à la
manière de calculer les toutes premières valeurs et comme il existe de nombreuses règles
empiriques, c'est la raison pour laquelle vous aurez peu de chance de trouver deux logiciels
qui donnent le même résultat pour ce modèle.
Nous pouvons obtenir si besoin est (très utile dans certains cas) que les coefficients du modèle
ou même un seul coefficient:
ou obtenir que les valeurs lissées des différentes composantes du modèle tel que calculé dans
le cours théorique:
Nous pouvons utiliser la commande native predict( ) pour faire des projections. Ici nous
allons faire un exemple aberrant allant au-delà d'une période:
Donc le meilleur modèle est un modèle ETS(A,A,N) ce qui signifierait que c'est un modèle
additif de Holt-Winters avec erreurs additives. La suite donne:
Donc ici nous avons un nombre d'informations très intéressantes avec les indicateurs d'erreurs
de mesures que nous avons vu dans le cours théorique. Graphiquement le meilleur modèle
donne:
À prendre évidemment avec des pincettes et toujours avec le même danger de lecture aberrant
que l'on retrouve dans toute la littérature spécialisée et que nous avons déjà mentionné dans le
cours théorique.
Nous allons ici calculer les coefficients d'autocorrélation avec les mêmes données que dans le
cours théorique après avoir donné la définition (intuitive) de l'autocorrélation.
Nous retrouvons donc exactement les mêmes valeurs que celles calculées à la main et que
dans le cours Minitab!
Nous allons ici calculer les coefficients d'autocorrélation avec les mêmes données que dans le
cours théorique après avoir donné la définition (intuitive) de l'autocorrélation.
Nous retrouvons donc exactement les mêmes valeurs que celles calculées à la main et que
dans le cours Minitab!
Donc selon R le meilleur modèle semble être un modèle AR(2). Graphiquement cela donne:
Les coefficients sont quand même significativement différents (environ 10%). Au niveau
graphique, cela nous donne quelque chose par contre de très semblable:
avec la command native ts.diag( ) nous pouvons avoir quelques informations équivalents sous
forme graphique et d'autres qui sont complémentaires:
Nous retrouvons donc bien la valeur de alpha vue plus haut (au signe près). Le graphique
associé donne:
et le graphique associé:
et le graphique associé:
Avec le package zoo et la commande aggregate par mois (as.yearmon), nous regroupons ces
données par mois:
Bon une fois ceci fait, nous chargeons le package forecast et utilisons la commande
auto.arima( ) qui va chercher le meilleur modèle arima possible (la recherche peut suivant la
Ensuite faisons un plot de ce que nous avons vu dans le cours théorique (sachant qu'il y a
beaucoup de choses que nous n'avons pas étudié non plus et que nous ne présenterons donc
pas ici):
Nous commençons avec la modèle AR(2) que nous avons obtenu plus haut). Comme la
commande native arima.sim( ) ignore l'interception, nous devons l'ajouter systématiquement:
et le graphique associé (qui donc au même titre que les valeurs sera différent à chaque re-
exécution des commandes ci-dessus):
et ainsi de suite... avec en fonction des modèles de subtilités pour rapatrier le constante donc
n'oubliez pas d'utiliser la commande str(data.arima) pour voir le contenu da la liste.
Nous partons donc des mêmes données que dans le cours théorique de 1267 lignes:
et le graphique associé:
Attention à ne pas oublier de contrôler avant les calculs que les résidus satisfont les conditions
d'application du modèle GARCH (voir le cours théorique)!
Alors avec cette première méthode que ce soit au niveau des paramètres ou du graphique cela
est très différent de toutes les méthodes théoriques que nous avons déterminées. Donc nous
laisserons de côté cette méthode à défaut de perdre du temps à décortiquer l'algorithme R
utilisé.
...
Alors que dire du résultat... D'abord la manière que R déterminer la moyenne empirique
nécessiterai de décortiquer le script car ce n'est pas la moyenne arithmétique, ni géométrique
ni la médiane...
Sinon pour les trois autres paramètres nous sommes à peu près bon et nous retrouvons les
mêmes résultats que dans le cours théorique que si nous considérons le t par rapport à la
moyenne empirique (alors que MATLAB le fait par rendement successifs).
En ce qui concerne le graphique, là aussi c'est très différent de calculs faits à la main dans le
cours théorique et aussi très différent de ce que nous avons obtenons avec MATLAB.
La transformée de Fourier Rapide est un outil très important dans l'étude des séries
chronologiques pour authentifier s'il y a des harmoniques périodiques dans une série de
mesures et la fréquence des ces dernières.
Je n'ai malheureusement pas trouvé de fonction ou de package proposant une FFT avec
exactement le même résultat que celui fourni par Matlab (qui reste souvent le logiciel étalon
pour ce qui est relatif au traitement du signal).
Donc à défaut voici une série de commandes (qui peuvent bien évidemment être mis sous la
forme d'une fonction simple) permettant d'obtenir exactement le même résultat que Matlab et
faisant appel à la fonctio native fft( ) de R (merci à un internaute anonyme!):
et voilà :-)
Certaines personnes préfèrent cependant les périodogrammes (dont je ne suis pas un fan):
Les données de Panel sont simplement un groupe des données longitudinales mesurées donc
pendant un certain temps. L'exemple le plus typique en entreprise étant la mesure des ventes
pendant plusieurs années d'articles (produits différents).
Voyons comment analyser ce type de données avec R. Nous chargeons d'abord le package
pscl qui contient un jeu de données bien adapté à ce type d'exemple (% de participaton aux
votations présidentielles américaines depuis 1932 par État et par éléction):
Ceci état pour poser le contexte. Ce qui nous intéresse est de comparer le tout. Donc déjà nous
pouvons superposer toutes les données:
Ce qui donne:
Nous savons que rien ne nous interdit de faire une ANOVA par État en passant par une
régression car nous avons étudié cela en cours.
La valeurde Fisher est significative ce qui indique que la régression par États peut être
considérée comme ayant tous ses coefficients significativement non nuls et donc que ce type
d'analyse à un intérêt certain.
Finance Quantitative
Nous allons d'abord commencer par la finance qualitative (simpler is better...). C'est-à-dire
l'affichage des valeurs d'indices ou d'actifs via de simples graphiques.
Nous utiliserons ici exclusivement le système de données de Yahoo (si vous souhaitez utiliser
le service Bloomberg voyez du côté de http://docenti.luiss.it/fasano/bloomr/). Si vous
cherchez à récupérer les données pour un indice ou une entreprise particulière dont vous ne
connaissez pas le code, vous pouvez utiliser la page suivante:
http://finance.yahoo.com/lookup
Pour commencer, nous allons utiliser le package ggplot2 et importer des données de Yahoo
en utilisant la commande geom_ribbon( ) pour faire une bande de fluctuation des cours par
rapport à la moyenne arithmétique:
Ce qui donne:
...et puis comme la finance c'est un peu jouer au casino (pile ou face) sur les marchés financier
un bon petit jeu est souvent ludique (et au bon passage une application ludique de la
commande sapply( )):
Nous détaillerons chacune des options avec un cas concret ultérieurement pour comparer les
résultats au cours théorique. Cela est important car le package GUIDE contient le minimum
qu'est censé connaître tout bon manager ou employé au bénéfice d'un MBA.
Attention!!! Lorsque vous changez les paramètres dans une boîte de dialogue, il faut valider
par la touche ENTER pour que les calculs prennent effet avec les nouveaux paramètres.
Nous allons pouvoir vérifier le bon fonctionnement de ce package avec par exemple un petit
calcul simple pour lequel nous ne trouvons aucune fonction dans MATLAB. Il s'agit du calcul
du prix future sur taux de change que nous avons calculé dans le cours théorique:
Nous allons dans l'exemple qui va suivre comparer au résultat obtenu avec le solveur de
Microsoft Excel et la méthode basée sur le norme de Frobenius mais en utilisant la fonction
nearPD( ) du package Matrix de R:
Bon je ne veux pas débattre et m'étendre sur ce modèle empirique mais juste présenter
brièvement sa mise en application car:
Bon étant donné ma passion dévorante pour ce modèle nous allons nous en tenir au minimum
syndical.
http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html
Pour télécharger les valeurs SMB (Small Minus Big: Facteur de taille) et HML (High Minus
Low: Facteur de prix) pour les U.S.:
S'y trouve alors un fichier avec les facteurs qui nous intéressent qui une fois un peu remané et
auquel on rajoute dans la dernière colonne notre investissement d'intérêt (il s'agit d'un Fond
comme un autre):
Après libre à vous d'utiliser le modèle de la même manière que pour la CAPM (voir page
1162). Bref la même chose peut être fait avec n'importe quel tableur d'entrée de gamme.
Dans le cours théorique de finance quantitative nous avons assez longuement étudié le
mouvement brownien. Nous avons cependant aussi parlé des processus de Lévy pour les
distributions à queues épaisses. Voyons-en un exemple en commençant d'abord avec le
processus gaussien:
et ensuite avec le cas classique du processus de Lévy basée sur ne distribution de Student
(remarquez le fait que parfois vous obtiendrez un grand saut):
et nous faisons quelques tests dans l'ordre de difficulté croissant au niveau du bagage
mathématique nécessaire que nous avons étudié dans le cours théorique d'économie
quantitative:
Nous pouvons obtenir l'équivalent avec la commande as.xts( ) du package xts associée à la
commande get.hist.quote( ) du package tseries:
et avec la commande getOptionChain( ) nous avons les informations concernant les options
(dérivés) du sous-jacent choisi (du moins quand certaines des informations ne sont pas
manquantes...):
http://finance.yahoo.com/q/op?s=GOOG&ql=1
Empiriquement j'ai mis un time de "1 minute" soit "60 secondes" qui affiche régulièrement
l'avancement dans la fenêtre de commande:
et:
p
log t log pt log pt 1 diff log i
pt 1
Mais sur une plus courte période (en-dessous d'un certain nombre d'années il y a un
décrochage et c'est assez drôle à voir), nous avons:
Ce qui donne:
Ce qui donne:
Si jamais voici un aperçu du fichier *.csv utilisé qui provient directement de Yahoo (s'inspirer
des exemples vus plus hauts pour importer des données de Yahoo!):
Nous ne voyons pas bien les prix OHLC (Open High Low Close). En restreignant la plage de
dates c'est mieux:
Si on veut enlever la partie concernant les volumes il suffit d'utiliser le paramètre TA:
ou encore:
ou encore:
Nous avoir que févirer 2011, nous écririons reChart(subset='201102'), pour avoir l'année
2010 jusqu'à juin 2011 nous écririons reChart(subset='2010/201106'), les dix premières
semaines de l'année choisie reChart(subset='first 10 weeks') ou les dix derniers mois
reChart(subset='last 10 months') ou le dernier mois zoomChart('last 1 months')
Si nous comparons le résultat avec le code VBA que nous avions écrit dans MS Excel, nous
voyons que nous avons pour les mêmes valeurs d'entrées la même valeur de sortie.
Voyons maintenant que le prix d'un Call pour une volalitlité donnée est bien une fonction
décroissante du strike (ce qui est assez logique):
Ce qui donne:
Et nous savons aussi que le prix d'un Call est aussi fonction croissant du temps à maturité:
Diagramme Payoff
Nous avons étudié dans le cours théorique les diagrammes payoff des Call et Put ainsi de de
plusieurs stratégies de portefeuilles d'options. Nous allons voir ici que ce n'est pas forcément
trivial à écrire.
C'est le diagramme payoff le plus simple que nous puissions imaginer et la syntaxe n'est pas
forcément triviale:
Ce qui donne:
Smile de volatilité
Rappatrions les données des options Microsoft négociables sur le marché de gré à gré (1
option est pour 100 sous-jacents normalement dans les exports Yahoo/Google) dont la date
d'expiration est avant le 30 Juillet 2014 et dans l'ordre décroissant des dates d'expiration
(comme le lecteur pourra le voir ci-dessous il n'y pas beaucoup d'options qui vont à une date
ultérieure au 19 Juillet 2014 sachant que ce script a été écrit le 1er Juillet 2014... donc il y peu
d'agents prenant le risque de proposer des options au-delà de 19 jours pour les sous-jacents
Microsoft):
Alors que la colonne Strike donne le... strike du sous-jacent (logique...), la colonne Last donne
le dernier prix de trade du Call. Nous nous baserons sur cette valeur pour déterminer la
volatilité implicite mais nous pourrions tout aussi bien prendre la valeur Bid (qui est pour
rappel le prix auquel les acheteurs esssaient d'acheter l'option) si nous sommes acheteurs de
Call.
La première colonne contient les dates d'expiration des Call (ou des Puts si vous descendez
plus bas dans la liste...) mais voici comment la décoder en détails:
Nous voyons que pour une même date d'émission des agents proposent des options d'achats
avec des strikes donnés variés (ce qui est normal puisque chacun fait des projections à sa
manière).
Nous allons pour poursuivre extraire la colonne des strikes, la convertir en valeurs
numériques, compter les nombres d'options disponibles dans cette extraction et déterminer et
extraire les noms de lignes pour ensuite les traiter afin d'obtenir les dates qui sont codées:
Ensuite nous calculons la durée jusqu'à maturité entre le jour de rédaction (exécution) du
script et la date d'expiration de chaque option (je n'ai pas trouvé autre chose qu'une boucle...):
Ensuite puisque par définition le smile de volatilité est pour des options de maturité égales,
nous allons filtrer sur un choix particulier arbitraire pour avoir que des options de même
maturité:
Il ne nous reste plus qu'à faire le graphique avec une interpolation de type spline:
Nous voyons que la smile (la volatilité implicite) est proche d'une minima pour la valeur du
cours du jour du sous-jacent. L'interprétation de cette courbe est très intéressante mais si nous
ne possédons pas toutes les informations du marché, nous pouvons a posteriori conclure que
le sous-jacent a une plus forte probabilité d'être à l'avenir au-dessus du cours du jour qu'en
dessous (c'est quand même des sous-jacents Microsoft....).
Surface de volatilité
L'extension naturelle du smile de volatilité est la surface de volatilité. Le code est quasiment
identique à avant à la différence que nous ne filtrons pas pour une seule date de maturité
donnée et aussi à la différence qu'il s'agit d'un plot 3D plutôt que 2D.
Ce qui donne:
Après ne reste plus qu'à faire un plot mais vu la faible quantité de données il ne va pas être
extra...:
Nous voyons donc bien ci-dessus la forme en "U" du smile de volatilité pour chaque maturité.
Bon cela doit être possible de faire plus joli... mais après quand on réfléchit à l'utilité de ce
type de graphiques 3D....
Volatilité implicité
Considérons maintenant une options Call qui le 13 Mai 2005 était évaluée à 0.0004 pour un
prix d'exercice K 23 , avec un spot S0 22.66 (visible sur le graphique ci-dessus en haut à
gauche). La date d'expiration était au 3 Juin 2005, ce qui correspond à 15 jours, dès lors
rapporté à une fraction d'année T 15 et le taux sans risque du marché était de
r 0.02074 . Dès lors, le prix de l'option était en utilisant un autre package que RQuantLib
(comme cela vous pouvez voir qu'il y en a plusieurs qui font la même chose):
La différence entre le prix du Call sur le marché est selon le modèle de Black & Scholes-
Merton est très différent. Voyons donc la volatilité implicite du marché:
Et ensuite, nous lançons un code de backtesting de stratégie (achat lorsque la valeur un peu
supérieure au SMA20):
Arbre binomial
Nous allons vérifier ici que nous obtenons bien avec R les calculs effectués à la main dans le
cours théorique concernant l'arbre binomial de Cox-Ross-Rubinstein.
Nous allons ici pouvoir appliquer ce nous avons étudié dans le cours théorique et comparer
les résultats – donc avec les mêmes données brutes – par rapport à l'approche approximative
que nous avions mis en place avec le solveur du tableur Microsoft Excel.
Ensuite, nous définissons combien de type d'actions nous voulons dans le portefeuille et
prenons empirique les 6 qui ont la plus grande valeur nominale et calculons la valeur globale
du portefeuille:
Nous construisons ensuite empiriquement un portefeuille avec les actions à poids égal
(remarquez bien la propriété type="equal"):
Nous pourrions également tout aussi bêtement choisir des poids proportionnels à la valeur de
chaque action dans le portefeuille en mettant la propriété type à linear:
Pour ceux qui aiment bien les graphes nous pouvons sortir un plot du résumé ci-dessus:
Et pour ceux qui aiment avoir le roic (Return On Invested Capital) sous forme de graph
tornado...:
Et concentrons-nous sur la ligne Technology. La colonne weigth est triviale. Retrouvons alors
la valeur de la colonne contrib. Le secteur Industrials est composé des trois actions suivantes:
Ensuite, nous regardons aussi l'exposition des différentes stratégies (on laisse le soin à
l'étudiant de vérifier que les calculs sont justes):
Considérons un gestionnaire de portefeuille d'actions qui utilise le S&P 500 comme indice de
référence. Un mois donné, il surperforme le S&P de 3% (le "rendement actif" pour rappel qui
est la différence entre le rendement du portefeuille et l'indice). Une partie de cette
performance s'explique par le fait qu'il a alloué plus poids du portefeuille à certains secteurs
qui ont bien performé. Appelez cet effet "l'effet d'allocation". Une partie de sa surperformance
s'explique par le fait que certains des stocks sélectionnés dans un secteur fait mieux que le
même secteur dans l'indice. Appelez ce l'effet "l'effet de sélection". Le résidu (la différence de
performance entre le rendement actif et les deux effets) peut alors être attribué à une
interaction entre répartition et la sélection - l'effet d'interaction.
Le modèle de Brisnon fournit des définitions mathématiques Brinson et des méthodes pour
calculer ces différents éléments.
wiB est le poids de chaque actif i dans l'indice (la somme étant égale à 1)
wiP est le poids de chaque actif j dans l'indice (la somme étant égale à 1)
W jB est le poids de chaque catégorie j dans l'indice:
W jB wiB , i j (la somme des W jB étant égale à 1)
i
Comme les poids dees catégories du portefeuille dynamiques sont généralement différents de
l'indice (portefeuille benchmark pour rappel), l'allocation joue un rôle dans le rendement actif
comme nous l'avons déjé mentionné. La même règle au niveau de la sélection des actifs dans
une catégorie.
N
Rallocation W jP W jB R Bj
j 1
où la multiplication par R Bj se justifie par le fait qu'indirectement nous supposons que dans
l'allocation que c'est uniquement le poids des catégories qui joue un rôle et non la sélection
des actifs dans chacune des catégories (donc le poids des actifs pour cahque catégorie est
supposé alors comme étant le même pour chaque portefeuille), raison pour laquelle pour
chacun des portefeuilles et pour une catégorie donnée, le rendement est alors supposé le
même et le rendement de référence est celui de l'indice.
N
Rsélection R Pj R Bj W jB
j 1
et là l'idée est similiare à celle du dessus à la différence que ce qui nous intéresse c'est pour un
poids de catégorie (en prenant celui de l'indice comme référence) d'analyser comme la
sélection des actifs à l'intérieur de celle-ci influence le rendement et donc in extenso c'est le
rendement de la catégorie qui va différer entre les deux portefeuilles.
Le rendement d'interaction se définit alors comme étant la part inexpliquée tel que:
parfois noté:
Faison un exemple:
Ensuite, nous utilisons la fonction brinson( ) pour obtenir un résumé de ce dont nous avons
parlé:
Modèle de Markowitz
D'abord nous chargeons le package fPortfolio et nous construisons le même petit portefeuille
de 3 actifs comme nous l'avons fait avec Microsoft Excel et ce exactement avec les mêmes
données:
Ensuite, nous allons spécifier quelques contraintes à notre portefeuille avec la commande
portfolioSpec( ). Voyons ce qu'il est possible de spécifier:
Une fois ces premières contraintes définies, nous allons spécifier déjà pour la suite si nous
voulons une stratégie Long (LongOnly), Short (Short) ou mixe (en spécifiant les vecteurs
comme ci-dessous):
Jusque là tout est toujours conforme à ce que nous avons fait dans le cours théorique! Nous
poursuivons donc!
Dans le cours théorique nous nous étions fixés un rendement précis à atteindre qui était de
22%. Pour trouver les poids actifs permettant d'atteindre cet objectif, nous utiliserons la
commande efficientPortfolio( ) du package fPortfolio:
Comme nous pouvons le voir, les résultats des poids est proche des calculs effectués à la main
dans le cours théorique.
Là encore nous allons vérifier si nous obtenons les mêmes poids que les calculs effectués dans
le cours théorique. Nous utilisons la commande minvariancePortfolio ( ) du package
fPortfolio:
Nous retrouvons donc exactement les mêmes poids que ceux calculés dans le cours théorique!
Frontière efficiente
Nous allons voir maintenant qu'il nous est possible de tracer la frontière efficiente avec la
commande frontierPlot( ):
En fixant le taux sans risque, nous pouvons pour tracer la courbe du ratio de Sharpe comme
nous l'avons vu dans le cours théorique avec la commande sharpeRatioLines( ):
Nous avons mis en évidence en gris l'ordonnée 0.22 avec l'axe en 0 de l'actif sans risque que
nous nous étions fixé dans le cours théorique (et là nous voyons que tout est bon puisque la
tangente passe bien par le point de coordonnée (0,0.22)).
Nous pouvons avoir un autre visuel résumant dans les grandes lignes les propriétés de notre
portefeuille en utilisant la commande tailoredFrontierPlot( ):
Avant de commencer il savoir que ce package ne peut manipuler que des séries temporelles de
type xts même si le contraire est écrit...
Nous nous concentrerons ici que sur des fonctionnalités du package dont les concepts sous-
jacents ont été démontrés dans le cours théorique et sur certaines fonctionnalités ne
nécessitant aucun package matématique particulier mais qui sont visuellement élégantes.
Et nous pouvons faire un plot des retours vis-à-vis d'un indice du marché:
Nous obtenons bien le même résultat qu'avec MATLAB. Le graphique associé nous donne:
Pour ceux qui veulent comprendre ce graph, voici le code source de la fonction (c'est
empirique et bête):
Commençons par l'indicateur un peu hors catégorie (car indépendant de toute probabilité dans
ce package) qui est l'écart-type seul avec un historique de 10 jours:
à 20 jours:
et ainsi de suite... Voyons maintenant la même chose avec les 3 VaR à un niveau de 95% que
nous avons vues dans le cours théorique avec aussi un historique à 10 et 20 jours:
et à 20 jours:
Pour ce faire, nous allons utiliser la commande ES( ) avec des données provenant aussi d'une
loi Normale:
MEDAF ou CAPM
Les données utilisées ici seront les mêmes que celles utilisées dans le cours théorique lors de
notre étude du modèle de diversification efficient de Sharpe.
Nous allons voir comment retrouver les paramètres du MEDAF comme beta, la prime de
risque… que nous avions calculé à la main.
D'abord calculons le rendement moyen de chacun des titres avec des commandes classiques
de R:
Donc jusque là les résultats sont parfaitement conformes à ce que nous avions calculé à la
main dans le cours théorique.
Comme nous l'avons démontré dans le cours c'est de l'algèbre élémentaire et donc nous
n'avions pas fait d'applications numériques. Faisons donc cela ici même exceptionnellement.
Rappelons que:
E Rm R f
Nous avons alors dans le cas présent pour le premier exemple du CAPM.RiskPremium( ):
Ensuite comme la rentabilité espérée était aussi triviale dans le cours théorique, nous n'avions
pas fait d'application numérique. Donc allons-y:
0.402 0.05
E RA 0.2178 0.05 0.87904051
0.03917
et donc nous n'obtenons pas la même chose. Cela semblerait à priori être une erreur du
package. Effectivement dans le code source, nous avons (version 1.1.0):
Donc il y a bien une erreur car à ma connaissance on ne doit pas multiplier par la moyenne
mais par l'écart-type!
Nous pouvons avoir des indicateurs de rations de performances avec la ration de Share
utilisant la commande SharpeRatio( ), l'alpha de Jensen avec la commande
CAPM.jensenAlpha( ) et enfin le ration de Treynor avec la commande TreynorRatio( ):
Bref rien d'extraordinaire ici... la finance est avec la statistique le domaine expert dans la
création à l'infini d'indicateurs empiriques...
Il existe aussi des fonctions pour calculer le résidu epsilon de la droite de régression et alpha
qui est l'ordonnée à l'origine en utilisant les commandes CAPM.alpha( ) et CAPM.epsilon( )
toujours du même package:
Ra f Rb Rb
Nous n'avions pas calculé dans le cours théorique l'epsilon pour la simple raison que c'est une
abstraction... L'auteur du package a fait le choix empirique et tout à fait discutable de
déterminer celui sur la base du calcul suivant:
n n
1 Ra 1 1 Rb 1
i 1 i 1
Donc il s'agit de calculer la différence entre les rendements composé des deux portefeuilles.
Choix tout à fait discutable...
i i2 m2 2
Risque total Risque systématique 2 Risque spcécifique 2
Pour laquelle nous n'avons pas fait d'application numérique pratique (état donné les
approximations implicites pour écrire cette relation qui par voie de conséquence est alors un
peu trop scolaire).
Value at Risk
Le but ici va être de vérifier si les calculs coïncident avec ce que nous avons vu dans le cours
théorique lors de notre étude des nombreux modèles de VaR.
Faisons maintenant un plot (malheureusement l'abscisse est imposée dans le code source du
package mais on comme on a accès au code source...) en utilisant la commande
chart.VaRSensitivity( ) et vérifions la correspondance avec les calculs faits précédemment
(raison pour laquelle j'ai rajouté une ligne verticale):
Commençons:
Ensuite:
1. La première ligne signifie que les éléments (séries chronologiques) 1 (A) et 2 (B) ont
été fusionnés et comme ils sont des singletons, le signe est négatif.
2. À la deuxième étape (deuxième ligne), nous avons les éléments (séries
chronologiques) 4 (D) et 5 (E) ont été fusionnés et comme ils sont des singletons, le
signe est négatif.
3. À la troisième ligne, comme les valeurs sont positives, cela signifie que nous avons
associé les deux premiers associations (donc celles des deux premières étapes)
ensemble.
4. À la quatrième linge, nous associons le 6ème (F) élément (6ème série chronologie) avec
l'association de la ligne 3 de cette même matrice {1,2}.
5. et ainsi de suite...
Nous poursuivons:
Ce qui donne à la première itération (puisque nous avons mis une pause avec la commande
readline( )):
à la deuxième itération:
à la troisième itération:
à la quatrième itération:
Pour la suite enlevez le readline( ) sinon quoi cela va poser des problèmes!!!
Maintenant, nous posons les séries chronologiques en face des feuilles ainsi que l'axe des
légendes:
plot(x[,r$order[i]],axes=FALSE,xlab="",ylab="",main="",type="l",col="navy",
lwd=2)
box()
}
par(op)
Actuariat
Le but de ce chapitre va être de vérifier les résultats renvoyés par le package
lifecontingencies correspond aux résultats obtenus suites aux quelques démonstrations
mathématiques triviales étudiées dans le cours théorique d'introduction élémentaire à la
finance actuarielle.
http://www.mortality.org
Nous chargeons donc d'abord le package lifecontingencies et le fichier *.csv qui nous
intéresse:
La colonne ex qui donne l'expérance de vie à la naissance est déjà plus intéressante à vérifier
pour voir si elle correspond aux notions vues dans le cours théorique. Pour obtenir l'espérance
de vie à la naissance relativement à cette table, nous savons qu'il suffit de faire la somme des
probabilités de survie de la fin jusqu'au temps d'intérêt divisé par la probabilité de survie
jusqu'à ce même tempd d'intérêt. Comme je n'ai pas trouvé dans R une commande unique
pour vérifier cela à la main nous pouvons dans le cas présent nous reporter à la feuille MS
Excel suivante:
Nous pouvons ensuite calculer la probabilité de la durée encore à vivre, c'est-à-dire dans le
cas de l'exemple ci-dessous la probabilité que si un individu est encore vivant à 20 ans qu'il
atteigne les 20+30=50 ans en utilisant la commande pxt( ):
Le package intègre plusieurs techniques empiriques d'interpolation pour des âges non entiers:
La probabilité de décéder pendant la même période sera donnée par la commande qxt( ):
Soit explicitement:
Il est aussi possible de travailler sur plusieurs têtes de plusieurs tables différentes de vie. Par
exemple si dans notre cas nous souhaitons savoir quelle est la probabilité jointe (pour rappel
c'est donc par définition le produit des probabilités) de survie que deux individus d'âge de
départ identiques ou différentes ont de survivre pendant la même période de temps, nous
utilisons alors la commande pxyt( ) qui correspond à h pxy :
Ou pour avoir la probabilité qu'une des personnes soit encore en vie (2 événements
incompatibles non nécessairement disjoints) h pxy :
Nous pouvons aussi contrôler que les calculs de capitals différés dans le cours théorique de
techniques de gestion sont accessibles avec ce même package lifecontingencies.
Calculons la valeur probable d'un capital qui pour rappel est donné par:
Sxm
m
Ex f em f em m px
Sx
Les calculs d'annuites sont alors probablement in extenso faciles à obtenir. Rappelons la
relation de la rente viagère temporaire praenumerando démontrée dans le cours théorique:
Bon je préfère utiliser le tableur Microsoft Excel pour cela mais suite à la demande d'un
participant à une de mes formations voici comment calculer le même VAN que dans le cours
théorique et qu'avec un tableur en utilisant la fonction presentValue( ):
Nous obtenons donc bien le même résultat que dans le cours théorique et qu'avec un tableur.
Où nous retrouvons bien la valeur de NFV (Net Future Value) vue calculée dans le cours
théorique:
Le package queuing a un commande résumant toutes les paramètres d'une file d'attente
M/M/1 tel que nous l'avons démontré dans le cours théorique en utilisant les commandes
QueueingModel( ) et NewInput.MM1( ):
Après nous pouvons faire cela pour de nombreuses autres files d'attentes avec ce package
mais nous nous limiterons poru l'instant à la M/M/1 car c'est la seule dont nous avons
démontré les propriétés mathématiques en détails.
Nous pouvons extraire les valeurs indépendamment les unes des autres (si nous avons besoin
de faire des graphiques par exemple):
Dans une entreprise, nous avons dénombré aux heures de pointe 200 appels d'une
duréemoyenne de 6 minutes à l'heure (temps de service moyen). Quelle est la probabilité de
saturation avec 20 opérateurs (sachant que la durée de service suit une loi exponentielle et la
distribution des arrivées une loi de Poisson)?
Quelle est la probabilité d'être mis en attente si nous prenons un nombre N de 7 serveurs s'il y
a un taux d'arrivée de 1 appel par minute et une durée moyenne de service de 5 minutes.
Statistiques bayésiennes
Le but ici n'est pas de redévelopper la discussion relativement au débat qui concerne les
statisticients fréquentistes contre les statistiques bayésiens mais uniquement de vérifier que
nous retrouvons bien les quelques rares calculs bayésiens fait à la main dans la cours
théorique ou ceux qui ont été simplement mentionnés avec suffisamment de détails
mathématiques.
Nous allons donc reprendre ici le même calcul dans le cours théorique basée sur l'étude des
cancers dans une école américaine a proximité de lignes haute tension.
Nous commençons par définir ce que nous savons et les hypothèses a priori:
Nous utilisons ensuite la commande pdisc( ) du package LearnBayes pour calculer les
probabilités a posteriori:
Nous retrouvons donc bien les probabilitées discrètes a posteriori calculées à la main dans le
cours théorique.
Nous avons démontré dans le cours théorique que si la distribution a priori est une distribution
Normale d'espérance 0 et de variance 02 et que la distribution du maximum de
vraisemblance est Normale, alors la fonction de distribution a posteriori aussi Normale
d'espérance a posteriori:
0 ny
02 2
1
1 n
2
0 2
et de variance a posteriori:
1
1 n
2 2
2
0
1
Supposons que nous souhaitons mesurer le quaotient intellectuel d'un group d'individus. Nous
avons une fort prémonation que a priori il 50% de probablité cumulée (médiane) que son QI
soit de 100 et qu'il y a 9/10 que celui-ci soit compris entre 80 et 120 (ces informations peuvent
aussi être basées sur des données historiques passées).
n 4, y 110
et que nous savons qu'au niveau mondial la variance vraie est donnée par 2 15 .
n4
y 110
2 15
0 100
02 12.159
0 ny
1
02 2 1 n
1 107.2442, 1 2 2 6.383469
1 n
2 0
0
2
L'idée ici va être d'estimer a priori la fonction de distribution de la proportion du cas discret
mais de façon continue. Comme il n'y a pas des milliers de fonctions de distributions dont le
support est compris entre ]0,1[ nous allons nous rabattre sur la seule que nous avons
rencontrée dans le cours théorique: la loi bêta.
Nous avons alors en faisant l'hypothèse que la médiane de la proportion de cancers est de 3%
et que le 66 centile est de 4% les paramètres de la fonction beta qui seront donnés par la
commande beta.select( ) du package LearnBayes:
Nous avons démontré dans le cours théorique que si la distribution a priori est une distribution
beta et que la distribution du maximum de vraisemblance est binomiale, alors la fonction de
distribution a posteriori est une distribution binomiale aussi avec les paramètres qui changent
un tout petit peu.
Ainsi, la probabilité a posteriori que la proportion soit plus grande ou égale à 3% est alors
donnée par:
Nous pesons un échantillon étalon de laboratoire sur une balance de précision d'écart-type
connu 3 et obtenons les masses suivantes:
Nous avons maintenant un nouvel échantillon et nous savons que a priori que sa masse devrait
être de moyenne 170 avec un écart-type aussi de 3 et distribué selon une loi Normale.
H 0 : 175, H1 : 175
H 0 P 175
H1 P 175
Donc l'hypothèse nulle a une cote a priori de 20 contre 1 d'être plus plausible que l'hypothèse
alternative. Et quid du a posteriori sachant que nous avons des mesures passées d'un
échantillon étalon?
n 10
y 176
2 9
0 170
02 9
0 ny
1
02 2 1 n
1 175.454545, 1 2 2 0.90453403
1 n
2 0
0
2
Ce qui nous donne les paramètres de la distribution a posteriori de la masse selon l'hypothèse
de Normalité.
Donc l'hypothèse nulle a une cote a posteriori de 1 contre 2 d'être plus plausible que
l'hypothèse alternative.
Bon et alors quoi? A priori l'hypothèse nulle est plus plausible mais a posteriori elle n'est le
pas! Eh bien nous allons devoir associer une probabilité à l'hypothèse nulle.
Suite à venir...
Statistiques spatiales
Une grande partie des statistiques spatiales sont relatives au fait de représenter graphiquement
des courbes de niveaux (isoclines) ou des densités surfaciques de probabilités ou simplement
des diagrammes à barre ou des diagrammes à secteurs representés sur des cartes
géographiques. Nous ne reviendrons pas sur tous ces cas déjà traités dans le chapitre sur les
Graphiques à partir de la page 463.
Pour cela nous générons des points aléatoirement qui nous serviront d'exemple:
Pour chaque point nous calculons la distance euclidienne aux autres points et nous traçons une
ligne pour chaque point avec son unique plus proche voisin (nn=nearest neighbour) ce qui
donne:
où le min(d[-i]) permet d'éliminer la distance nulle entre le point traité dans la boucle et lui-
même.
Ce qui donne:
Ce qui donne:
Nous pouvons aussi calculer "à la main" ou avec le package de statistiques spatiales spatstat et
sa fonction nndist( ) la distance moyenne:
Nous avons démontré dans le cours deque l'espérance de la distance du plus proche voisin
dans le cas d'une distribution aléatoire était donnée dans un disque par:
1
r E R
2
Or nous avons dans notre exemple 100 points dans un carré de côté unitaire. En première
approximation avec un disque inscrit au carré, nous avons alors:
100
127
R2
Dès lors:
1 1
r E R 0.044
2 2 127
Comme 0.044 n'est pas trop éloigné de 0.0500 la distribution peut être considérée
effectivement comme aléatoire.
Text Mining
Le texte mining est un sujet qui a toujours été important dans certains domaines académiques
mais depuis l'émergeance des réseaux sociaux et de l'analyse de news en temps réel, ces
techniques ont connues une véritable explosion.
Dans ce chapitre nous allons nous concentrer sur le techniques proposés par le par le package
tm (text mining) de R en se basant sur un corpus contenant un seul et unique fichier qui est le
pdf de mon e-book Éléments de Mathématiques appliquées (~4'900 pages A4) convertit en
texte pur (*.txt) ce qui le réduit à environ 2'800 pages A4.
Supprimer la ponctuation:
Il y a possibilité de supprimer aussi les mots utlra courants d'une langue donnée et qui
n'apportent rien à une analyse. Voyons d'abord quels sont ces mots dans le cas de la langue
française:
ou nous pouvons aussi enlever nos propres termes (en réalité il vaut mieux avoir un fichier
externe avec l'ensemble des termes considérés comme inutiles car en français il y en a un tel
nombre que la commande stopwords utilisée précédemment est quasiment inutile):
Arrivé à ce point nous ne pouvons pas encore lancer une analyse du corpus. Nous allons
devoir effecteur un autre nettoyage (voir ci-après):
Pour ne pas avoir à nettoyer cela manuellement (ce travail étant laborieux). Nous pouvons le
laisser au package SnowballC mais mixé avec la commande tm_map( ) du package tm (pas
besoin de spécifier la langue il s'en occuppe) cependant le résultat du travail est à ce jour
(2014-07-07) mauvais en français:
Une fois le nettoyage terminé, pour analyser le corpus nous utilisons la commande
DocumentTermMatrix( ) du package tm:
Voici par exemple le temps nécessaire pour analyser la version texte brute de mon livre (pour
que le lecteur ait une idée du temps de traitement de fichiers volumineux):
Ensuite, nous affichons un résumé de la variable créée pour voir le nombre de termes
analysés:
Ensuite, nous transformons la variable dtm en tant que matrice et faisons la somme des
colonnes afin d'avoir la fréquence d'apparition de chaque mot:
Ensuite, avec les commandes classiques natives de R nous pouvons trier et afficher les mots
les plus fréquents (remarquez les apostprophes qui n'apparaissent pas car le package tm ne
gère par les apostrophes ' mais veut une apostrophe stylisée ´):
Nous pouvons détecter les langues avec la commande textcat du package textcat( ) nous
avons avec un exemple simpliste:
Ou en utilisant notre texte de tout à l'heure (j'ai aussi mis le temps afin que le lecteur puisse se
faire une idée):
Bioinformatique/Biostatistiques
La bioninformatique (dont nous retrouvons aussi des Toolbox dans MATLAB) est aussi
excellent exemple d'application pratique des outils statistiques et de machine learning.
La technologie des microarray utilise les propriétés de l'hybridation des acides nucléiques et
utilise des molécules complémentaires attachées à une surface solide, nommées "probes" dans
le domaine, pour mesurer la quantité d'une acide nucléique spécifique (ou d'une chaîne
particulière) présent dans un échantillon souvent nommé "cible".
Les molécules se retrouvant dans les "probes" sont par la suite labelissées à l'aide de scanners
spécialisés et la quantité d'acides nucléiques dans chaque probe est donnée en termes
"d'intensité" (fluorescence de l'activié chimique) ce qui se traduit souvent sur les ordinateurs
en une image où chaque "probe" est un pixel ou ensemble constant de pixels ayant un couleur
donnée avec un certaine intensité (souvent du noir – rien - au blanc pour le maximum
d'intensité). La matrice (grille) de couleurs obtenue se lit parfois en ligne pour un type de gène
donné et en colonne des conditions d'analyse du gène d'intérêt.
Les industriels12 proposent une large palette de plateformes d'analyses différentes dont le
critère majeur est souvent le bruit de mesure (et en deuxième le type de cible d'intérêt:
génome humain, rat ou autre parfois sur mesure…) qui dans la pratique nécessitent des
répétitions d'analyses et manipulations statistiques des données extraites avant de pouvoir
conclure sur quoi que ce soit de manière un peu scientifique (donc pour un même échantillon
on obtient des résultats sensiblement différents à cause du bruit, du protocole expérimental et
de l'opérateur). Ces manipulations sont souvent nommées "preprocessing" (prétraitement)
dans le domaine.
12
En ce début de 21ème siècle Affymetrix serait le leader du marché des microarray.
Nous avons aussi besoin du package seqinr pour obtenir les séquences ADN et extraire des
informations des bases de données de protéines en lignes et aussi pour leur analyse
élémentaire:
Ces trois échangent des données des bases de données chaque nuit, donc en un point
quelconque dans le temps, elles contiennent des données presque identiques.
Chaque séquence dans la base de données NCBI séquence est mémorisée dans un registre
séparé, et se voit attribuer un identifiant unique qui peut être utilisé pour se référer à cet
enregistrement de la séquence. L'identifiant est connu comme une adhésion, et consiste en une
combinaison de chiffres et de lettres.
Par exemple, le virus de la dengue (grippe tropicale) provoque de la fièvre de la dengue, qui
est classée comme une maladie tropicale et désignée par l'un de quatre types de virus de la
dengue: DEN-1, DEN-2, DEN-3 et DEN-4. Les adhésions NCBI pour les séquences d'ADN
de la DEN-1, DEN-2, DEN-3 et DEN-4 sont respectivement NC_001477, NC_001474,
NC_001475 et NC_002640.
Note that because the NCBI Sequence Database, the EMBL Sequence Database, and DDBJ
exchange data every night, the DEN-1 (and DEN-2, DEN-3, DEN-4) Dengue virus sequence
will be present in all three databases, but it will have different accessions in each database, as
they each use their own numbering systems for referring to
Notez que parce que la base de données de séquences NCBI, EMBL et DDBJ synchronises
leurs tous les soirs, les séquences du virus DEN-1 (et DEN-2, DEN-3, DEN-4) sera présent
Nous allons voir maintenant comment extraire la séquence NC_001477 du NCBI par
exemple. Pour cela nous allons sur le siter Internet du NCBI et nous tapons le séquence qui
nous intéresse dans le champ de recherche:
Cliquez sur Create File, pour que ce fichier texte à l'extension *.FASTA se télécharge sur
votre ordinateur. Vous aurez alors:
Nous allons voir maintenant comment faire quelques petites traitement et analyses
élémentaires du fichier FASTA que nous avons téléchargé.
ou aussi:
Il se semble qu'il y ait a priori une anomalie avec le microarray h. Nous pouvons investiguer
cela un peu mieux avec un histogramme:
On voit qu'un des microarray a une distribution bimodale ce qui est souvent le signe d'un
défaut de conception.
Si nous décompressons un de ces fichiers, nous avons alors à chaque fois un fichier d'environ
1MB:
Nous pouvons avoir un petit résumé simplifié de notre jeu de fichier *.CEL:
Nous allons ici présenter la syntaxe pour résoudre et représenter les équations différentielles
que nous avons démontrées dans le cours théorique. À chaque équation différentielle sera
dédiée un exercice.
À une variante de notation près au niveau des variables par rapport au cours théorique nous
avons pour ce système d'ODE:
Ce qui donne:
À une variante de notation près au niveau des variables par rapport au cours théorique nous
avons pour ce système d'ODE:
Ce qui donne:
Nous tirons à pile ou face et regardons s'il y égalité entre les piles et les faces:
Souvent il n'y aura pas égalité (logique) mais nous savons cependant que la pièce n'est pas
truquée. Alors bon combien de fois avons nous la non égalité si nous répétons souvent ce test?
Donc 74 fois sur 1'000 l'hypothèse nulle comme quoi la pièce est non truquée est juste, le
reste du temps elle est fausse. Donc nous avons un risque de première espèce d'environ
93% de rejeter à tort l'hypothèse nulle.
Donc sur 1'000 expérience d'une pièce truquée, notre test à décidé 7 fois à tort que l'hypothèse
nulle était vraie à tort et 993 fois qu'elle fausse à raison. Nous avons donc un risque de
deuxième espèce de 0.007%
Dans tout language de script ou de programmation la notion de portée des variables est
importante à comprendre et dans certains cas très importante à utiliser.
Nous remarquons que la fonction est capable de lire la variable y qui est à l'extérieur de son
domaine de définition bien que cette même variable ne soit pas passée en tant qu'argument.
Ensuite, nous remarquons aussi que la fonction n'est pas capable de modifier la variable x qui
se trouve à l'extérieur de son domaine de définition puisque avant et après x vaut toujours 3.
Il faut voir maintenant comment nous pouvons écraser la variable x au-delà de la fonction
elle-même. Pour cela, la première méthode consiste à utiliser l'opérateur de superaffectation
<<-:
Et nous écrivons:
Ensuite, nous enregistrons notre fichier en tant que fichier script (pouvant contenir plus qu'une
fonction):
Nous un premier exemple avec la structure conditionnelle ifelse( ) qui est très souvent utilisé
pour opérer directement sur l'ensemble des attributs d'un vecteur:
Le ifelse peut aussi être utilisé pour autre chose que de faire des caluls bien évidemment (par
exemple renommer le contenu d'un vecteur-colonne) et ainsi de suite...
Voyons un petit cas classique qui consiste à faire un petit menu à choix exécutant certaines
commandes en fonction du choix de l'utilisateur:
Les structures itératives ont elles aussi une place très important dans R (comme dans tous les
langages de toute façon!). Considérons le cas où nous souhaiterions avec une identification
rapide des lignes paires ou impaires d'un data frame:
Ensuite, nous enregistrons notre fichier en tant que fichier script (pouvant contenir plus qu'une
fonction):
Ou encore l'instruction while( ) qui permet d'identifier par exemple à partir de quel moment
une certaine somme est atteinte:
Nous pouvons nous amuser à en faire une fonction avec un paramètre d'entrée:
Ce qui donnera:
R se base principalement sur la manipulation de vecteurs et matrices donc plutôt que de faire
des boucles parfois inutilement pensez vectoriel!!! Un exemple typique consistant à calculer
la somme des carrés des erreurs d'un vecteur:
Comme premier script, nous allons écrire une fonction qui permet de calculer au choix la
moyenne arithmétique, géométrique (ce qui nous évitera de faire appel au package psych qui
contient une commande geometric.mean( ) pour la moyenne géométrique), harmonique et la
médiane.
et nous écrivons notre script pour découvrir par la même occasion la structure conditionnelle
switch:
Ensuite, nous enregistrons notre fichier en tant que fichier script (pouvant contenir plus qu'une
fonction):
Remarque: Certes ce bout de code est discutable car s'il y a une valeur négative ou nulle dans
le vecteur, la moyenne géométrique va renvoyer une erreur (voir le cours théorique pour
traiter de ces cas là).
Nous écrivons (en même temps c'est une jolie application de la commande sapply( )):
et nous testons:
Programmer une fonction R avec deux paramètres en entrée est très similaire à n'importe quel
autre langage. Ce qui est drôle c'est de passer un vecteur à un des deux paramètres.
Considérons l'exemple physique avec le couple Terre/Soleil et l'apogée et à la périgée:
Ou une fonction qui peut être utile dans la pratique en plus d'être instructive (à personnaliser
selon les besoins) qui va authentifier dans un vecteur les positions d'une séquence continue de
k éléments identiques (dans la cas présent il s'agit de "1"):
Il s'agit ici probablement du cas le plus courant dans les complications pratiques du logiciel.
Écrire une fonction qui prend un fichier texte en entrée est identique à écrire une fonction qui
prend en entrée un scalaire ou un vecteur.
Pour cet exemple, nous allons prendre un fichier texte simple et générer un objet de type liste
dont chaque entrée correspondra à un mot se trouvant dans le fichier texte et nous associerons
à chacune de ses entrées, un vecteur contenant la position des endroits où ce même mot a été
trouvé.
La fonction dans un premier temps est la suivante (prenez le temps de lire car c'est instructif!):
et ainsi de suite avec toutes les techniques déjà étudiées pour les listes au début de cet e-book.
Un grand classique qui a son utilité pour être étendu à de nombreux autres domaines dans la
pratique:
Ou une façon plus amusante avec un indicateur d'avancement de calcul (très utilisé dans la
pratique) combinant les commandes winProgressBar( ) et Sys.sleep( ):
Signalons encore une autre approche qui montre la convergence du calcul mais utilisant cette
fois-ci le disque en entier (et pas que le premier quadrant):
Ce qui donne:
Dans les grosses applications, l'optimisation de scripts est un élément crucial (c'est un métier
en soit par ailleurs). Voyons un joli petit exemple avec R qui consiste à calculer les
différences (variations) éléments par éléments d'un vecteur par ordre d'apparition (certes il
existe un package qui contient déjà cela mais bon c'est instructif quand même!).
Maintenant faisons le code "au pire" qui consiste à générer un vecteur des différences avec
une simple boucle:
Comme nous pouvons le voir... c'est même plus lent qu'un tableur...
Mais nous pouvons faire encore mieux en utilisant les capacités vectorielles de R:
Statistiquement nous pouvons aller encore plus vite en utiisant le package natif compiler de R
et sa fonction cmpfun:
Ensuite il n'y a qu'à jouer avec en exécutant la fonction et ensuite en l'appelant avec les deux
dernières lignes du script visible ci-dessus.
Un grand nombre d'entreprises utilisent MS Excel pour faire des "statistiques". L'outil étant
très limité, une possibilité très intéressant est de faire des routines VBA qui exportent les
données en *.csv, laisser R faire les analyses et ensuite exporter les résultats ou graphiques
sous forme d'image pour les importer enfin dans MS Excel. Nous allons ici juste montrer le
principe qui se généralise ensuite à des cas beaucoup plus complexes. Pour plus
d'informations, je renvoie le lecteur à mon e-book sur le VBA dans MS Excel.
D'abord, nous considérons le script suivant RVBA.R enregistré dans le dossier C:/tmp/:
et ensuite dans un module VBA d'un des logiciels de la suite Microsoft Office, nous écrivons:
Lorsque l'on traite plus d'une dizaine de millions de lignes de données, travailler avec les
outils standard de R devient chose impossible et surtout si aucun package ne correspond à nos
besoin. Dès lors, soit on part sur du Matlab pour compiler un script en C, soit on écrit
directement du C dans R.
Pour voir cela avec un exemple simple et générique, nous allons reprendre l'exemple
générique donné initialement ici:
http://dirk.eddelbuettel.com/blog/2012/11/20/
http://cran.r-project.org/bin/windows/Rtools/
et nous pouvons exécuter notre petit script qui va donc estimer la valeur de Pi par Monte-
Carlo en utilisant le package Rcpp qui découle de l'installation des R Tools ainsi que le
package rbenchmark qu'il faudra avoir installé après coup et qui est spécialisé dans la
comparaison de scripts:
Le C# est un langage très demandé dans la finance dans certains pays. Nous allons donc voir
ici une méthode simple pour envoyer à R des valeurs et récupérer des résultats.
Ensuite ouvrez Visual Studio C# (ici la version 2010 dans un environnement XP x86) et créez
un Nouveau projet…:
Une fois le projet créé, cliquer sur Ajouter une référence dans menu contextuel du volet
droit sur le dossier Références:
Ensuite il est conseille dans certains environnements Microsoft Windows d'ajouter le chemin
de R dans la variable d'environnement PATH:
Une fois ceci fait, nous écrivons le code suivant pour faire un exemple d'appel de R:
using Microsoft.Win32;
using System;
using System.Linq;
using RDotNet;
using System.IO;
namespace Sample1
{
class Program
{
static void Main(string[] args)
{
using (RegistryKey registryKey =
Registry.LocalMachine.OpenSubKey(@"SOFTWARE\R-core\R"))
{
var envPath = Environment.GetEnvironmentVariable("PATH");
string rBinPath = (string)registryKey.GetValue("InstallPath");
string rVersion = (string)registryKey.GetValue("Current Version");
rBinPath = System.Environment.Is64BitProcess ? rBinPath + "\\bin\\x64"
:
rBinPath +
"\\bin\\i386";
Environment.SetEnvironmentVariable("PATH",
envPath + Path.PathSeparator + rBinPath);
}
using (REngine engine = REngine.CreateInstance("RDotNet"))
{
// Initializes settings.
engine.Initialize(); // After Executing this line its crashing.
}
}
}
Le but ici va être de découvrir des commandes qui sont similaires dans l'idée aux techniques
de debuggage de MATLAB:
Traceback
La commance traceback affiche la liste des fonctions appelées avant de lister enfin ce qui a
fait bugger le traitement.
Nous voyons que la commande traceback( ) est finalement a priori peu utile car renvoie des
informations trop vagues (du moins sur l'environnement MS Windows). Effectivement, si
nous changeons l'emplacement de l'erreur, nous avons...:
Donc la fonction h(y) apparaît aussi alors qu'elle ne génère aucun bug directement!
Debug
Bon nous remettons le code buggé dans l'état initial et nous allons voir maintenant la fonction
debug( ) dont l'utilité est aussi discutable. Effectivement voyons un exemple:
Comme nous pouvons le voir il est donc important de connaître (identifier) le nom des
fonctions bugées au fur et à mesure que l'on avance. Ainsi, dans l'exemple suivant, nous
décidons de débugger la fonction h(z) pendant le débuggage de g(y):
n: poursuivre le debuggage (sans que cela rente nécessairement dans les sous
fonctions!)
c: exécute l'ensemble restant de la fonction (ou de l'ensemble des fonctions) jusqu'à la
fin
Q: arrête le débuggage et l'exécution de la fonction
print(nom_variable): pour pouvoir à tout moment voir le contenu (la valeur) d'une
variable pendant un débuggage
Remarque: Il suffit de recharger le fichier *.R ou le package pour que le mode débuggage
qui a été activé en écrivant debug(nom_de_la_fonction) s'arrête
Browse
Il existe une autre manière beaucoup plus agérable de débugger qui consiste à écrire la
commande browse( ) à un endroit désiré du code source:
Comme vous pouvez le voir le debuggage commence donc à l'endroit où nous avons
positionné la function browse( ).
Trace
Il est possible lorsque les fonctions sont petites de les modifier "temporairement" pour faire
des tests à la volée pendant un debuggage. Pour voir cela considérons toujours notre petit
script:
Profilage
Le profilage de script est un sujet important pour ceux cherchant à trouver les goulets
d'étranglement de leurs codes. Pour l'exemple, nous utiliserons le même script qu'avant avec
les fonctions natives Rprof( ), invisible( ) dont le but est de ne pas imprimer les calculs de la
fonction testée, et summaryRprof( ).
Avec un tout petit script comme celui que nous avons, nous voyons que le profiler a de la
peine à profiler quoi que ce soit:
D'abord pour faire rapidement (et mal) un petite package R d'un script R, l'utilisatin du
package devtools sera obligatoire:
Une fois ce package installé, vérifiez votre dossier de travail par défaut et créez un dossier
pour un futur package avec la commande create( ) du package devtools:
Ce qui donnera:
Mais une erreur va empêcher l'installation. C'est parce qu'il fallait ajouter un fichier
NAMESPACE au préalable avant de compiler:
Nous pouvons ensuite utiliser notre fonction comme toute autre fonction:
Autre exemple:
Ou:
doplot = function(...){
m = as.numeric(tclvalue(mean))
s = as.numeric(tclvalue(sd))
n = as.numeric(tclvalue(n))
genplot(m,s,n)
}
base = tktoplevel()
tkwm.title(base,'Simulation Loi Normale')
mainfrm = tkframe(base)
mean = tclVar(10)
sd = tclVar(1)
n = tclVar(1000)
img = tkrplot(mainfrm,doplot)
scalefunc = function(...)tkrreplot(img)
sc_1 = tkscale(mainfrm,command=scalefunc,from=-10,to=10,showvalue=TRUE,
variable=mean,resolution=.01,orient='horiz')
m = tkframe(mainfrm)
tkpack(tklabel(m,text='Mean, SD',width=10),side='left')
tkpack(tkentry(m,width=5,textvariable=mean),side='left')
tkpack(tkentry(m,width=5,textvariable=sd),side='left')
tkpack(m,side='top')
tkpack(sc_1)
tkpack(sc_2)
tkpack(sc_3)
tkpack(mainfrm)
Ce qui donne:
http://rcom.univie.ac.at/download/current/statconnDCOM.latest.exe
Une fois le fichier télécharé vous lancez l'installation ce qui donnera sur un environnement
Microsoft Windows:
et Next...:
et Next:
et Next:
et on clique sur Finish. Une fois ceci fait, vous télécharger le logiciel RExcel sur le lien
suivant:
http://rcom.univie.ac.at/download.html
Ci-dessus vous cliquez sur I accept the agreement et ensuite sur Next:
Vous lisez très attentivement ce qui est écrit en rouge et cliquez sur Next:
Encore sur Next et apparaît le message suivant qu'il faut lire attentivement!:
Vous cliquez alors sur Oui et apparaît le second message suivant qu'il faut lire attentivement!:
Normalement vous aurez alors sur le bureau l'icône spéciale RExcel... qui aura été rajoutée
automatiquement:
Le minimum mimorum est de savoir convertir une plage de données Microsoft Excel en un
data frame R. Pour cela, après avoir saisi des données, on fait un clic droit et on choisit Put R
DataFrame:
Ensuite, nous pouvons vérifier que ce dataframe a bien été créé en faisant une petite analyse
élémentaire:
Il vient alors la boîte de dialogue ci-dessous dans laquelle il faut sélectionner en bleu les
variables voulues pour que RCmdr sache ce qu'il doit prendre comme variables en entrée et en
sortie:
Pour cela, il suffit d'utiliser d'appeler la function RApply( ) en passant en argument le nom de
la fonction R et ensuite les paramètres correspondants:
unz [UNZ]ip
[GET] [W]orking
getwd
[D]irectory
dir [DIR]ectory
[S]tring [PRINT]
sprintf
[F]ormatted
c [C]ombine values
regexpr [REG]ular [EXPR]ession Why "regular"? See regular sets, regular language
[DIAG]onal values of a
diag
matrix
col [COL]umn
lapply [L]ist [APPLY] Apply function to each element and return a list
MARGIN = 1 or rows [1] come before e.g., a 2 x 3 matrix has 2 rows and 3 columns (note:
2 in apply columns [2] row count is stated first)
cor [COR]relation
[AN]alysis [O]f
ancova
[COVA]riance
[M]ultivariate [AN]alysis
manova
[O]f [VA]riance
[T]ukey's [H]onestly
TukeyHSD
[S]ignificant [D]ifference
[H]ierarchical [CLUST]er
hclust
analysis
[C]lassical metric
cmdscale [M]ulti[D]imensional
[SCAL]ing
[PRIN]cipal [COMP]onents
princomp
analysis
[PR]incipal [COMP]onents
prcomp
analysis
resid [RESID]uals
[V]ariance-[COV]ariance
vcov
matrix
[B]ayesian [I]nformation
BIC
[C]riteria
[EVAL]uate an R
eval
expression
[READ] a file in [C]omma i.e., in each row of the data commas separate
read.csv
[S]eperated [V]alues format values for each variable
attr [ATTR]ibute
rev [REV]erse
prod [PROD]uct
var [VAR]iance
sd [S]tandard [D]eviation
intersect [INTERSECT]ion
[IM]aginary part of a
Im
number
[T]ranspose of a vector or
t
matrix
substr [SUBSTR]ing
[SUB]stitute identified
sub
pattern found in string
[G]lobal [SUB]stitute
gsub identified pattern found in
string
[N]umber of [CHAR]acters
nchar
in a string
[WIN]dows [METAFILE]
win.metafile
graphic
hist [HIST]ogram
[PLOT] colums of
matplot
[MAT]rices
[Q]uantile-[Q]uantile [P]lot
qqnorm
based on normal distribution
h argument in
[H]orizontal line
abline
v argument in
[V]ertical line
abline
[C]haracter [EX]tension or
cex as par [EX]pansion of plotting
objects
[C]haracter [EX]tension or
cex.sub as par
[EX]pansion of [SUB]title
[C]haracter [EX]tension or
cex.axis as
[EX]pansion of [AXIS]
par
annotation
[C]haracter [EX]tension or
cex.lab as par [EX]pansion X and Y
[LAB]els
xyplot [X] [Y] [PLOT] [X] for horizontal axis; [Y] for vertical axis
qq [Q]uantile-[Quantile] plot'
optim [OPTIM]isation
lm [L]inear [M]odel
[G]eneralised [L]inear
glm
[M]odel
[LO]cally [E]stimated
loess
[S]catterplot [S]moothing
[R]andom number
rexp
generation from
[EXP]onential distribution
[R]andom number
rgamma generation from [GAMMA]
distribution
[R]andom number
rpois generation from [POIS]on
distribution
[R]andom number
rweibull generation from
[WEIBULL] distribution
[R]andom number
rcauchy generation from [CAUCHY]
distribution
[R]andom number
rbeta generation from [BETA]
distribution
[R]andom number
rt generation from [t]
distribution
[R]andom number
rf generation from [F] F for Ronald [F]isher
distribution
[R]andom number
rchisq generation from [CHI]
[SQ]uare distribution
[R]andom number
rbinom generation from
[BINOM]ial distribution
[R]andom number
rgeom generation from
[EXP]onential distribution
[R]andom number
generation from
rhyper
[HYPER]geometric
distribution
[R]andom number
rlogis generation from [LOGIS]tic
distribution
[R]andom number
rnbinom generation from [N]egative
[BINOM]ial distribution
[R]andom number
runif generation from [UNIF]orm
distribution
[R]andom number
rwilcox generation from
[WILCOX]on distribution
geom_ in
[GEOM]etric object
ggplot2
stat_ in
[STAT]istical summary
ggplot2
coord_ in
[COORD]inate system
ggplot2
qplot in
[Q]uick [PLOT]
ggplot2
FUN as
[FUN]ction
argument
pos as
[POS]ition
argument
lib.loc in
[LIB]rary folder [LOC]ation
library
sep as
[SEP]erator character
argument
comment.char [COMMENT]
in read.table [CHAR]acter(s)
[I]nhibit [I]nterpretation or
I
[I]nsulate
T value [T]rue
F value [F]alse
[M]edian [A]bsolute
mad
[D]eviation
Show [STR]ucture of
ls.str
[L]i[S]ted objects
envir as
[ENVIR]onment
argument
q [Q]uit
ave [AVE]rage
min [MIN]imum
max [MAX]imum
[N]umber of [LEVELS] in a
nlevels
factor
[G]eneralised [L]east
gls
[S]quares
dwtest in
[D]urbin-[W]atson Test
lmtest
[S]tructural [E]quation
sem in sem
[M]odel
betareg in
[BETA] [REG]ression
betareg
[D]egrees of [F]reedom of
df.residual
the [R]esidual
deriv [DERIV]ative
[CHOL]eski [2=TO]
chol2inv
[INV]erse
[S]ingular [V]alue
svd
[D]ecomposition
[EIGEN]value or
eigen
[EIGEN]vector
[LOWER] [TRI]angle of a
lower.tri
matrix
[UPPER] [TRI]angle of a
upper.tri
matrix
[A]uto [C]orrelation or
acf
[C]ovariance [F]unction
[C]ross [C]orrelation or
ccf
[C]ovariance [F]unction
StatET as
Anyone know? Statistics Eclipse?
software
JGR as
[J]ava [G]UI for [R] pronounced "Jaguar" like the cat
software
ESS as
[E]macs [S]peaks [S]tatistics
software
Rterm as
[R] [TERM]inal
program
R CMD as
I think: [R] [C]o[m]man[D] prompt
program
repos as
[REPOS]itory locations
option
.site file
[SITE] specific file e.g., RProfile.site
extension
Frank [HARRELL]'s
Hmisc package package of [MISC]elaneous
functions
c in debug [C]ontinue
Q in debug [Q]uit
[PSYCH]ology related
psych package
functions
strucchange
[STRUC]tural [CHANGE]
package
relaimpo
[RELA]tive [IMPO]rtance
package
[C]ompanion to [A]pplied
car package Named after book by John Fox
[R]egression
df in
[D]ata [F]rame
write.foreign
R [F]requently [A]sked
R FAQ word
[Q]uestions
[D]e[V]ice [I]ndependent
DVI format
file format
[M]aximum [L]ikelihood
mle
[E]stimation
[S]ocial [N]etwork
sna package
[A]nalysis
[E]xponential [R]andom
ergm package
[G]raph [M]odels
Lorsqu'on décortique un script R dont on est pas l'auteur, souvent il nécessaire d'identifier le
type d'object que l'on manipule et parfois de la convertir en un autre type d'objects. Nous en
verrons quelques cas d'application concrets par la suite mais sinon voici au besoin la liste des
commandes utiles pour authentifier le typage et faire une conversion: