The French version of this document is no longer maintained: be sure to check the more up-to-date English version.
Spécificités de R
Exemples
R : documentation
ESS
R : quelques fonctions élémentaires
R : structures de données
Opérateurs
R : structures de contrôle
R : fonctions
Dans cette partie, nous faisons un rapide survol du logiciel : où se situe-t-il face à d'autres logiciels de calculs ? Que peut-on faire avec ? Quelles sont ses limitations ? Quelle est sa syntaxe ?
Il s'agit d'un logiciel de statisique : contrairement à d'autres logiciels de calcul numérique (Scilab, Octave),
http://www-rocq.inria.fr/scilab/ http://www.octave.org/
il y a déjà des fonctions pour effectuer des traitements statistiques non triviaux, qu'il soient classiques (régression, régression logistique, analyse de la variance, arbres de décision, analyse en composantes principales, etc.) ou plus modernes (réseaux de neurones, bootstrap, etc.).
Il s'agit d'un vrai langage de programmation, et pas un cliquodrome : on n'est pas limité par l'imagination des concepteurs du logiciel, on peut l'utiliser à n'importe quelle fin.
Il s'agit d'un langage interprété : l'avantage, c'est qu'on passe moins de temps à programmer avec, l'inconvénient, c'est que les calculs sont plus lents -- mais à moins de vouloir faire des calculs en temps réel sur les fichiers de la sécurité sociale, en tenant compte de toute la population française, c'est suffisemment rapide. Si on réellement besoin de rapidité, on peut se tourner vers SAS (commercial) ou DAP (libre).
http://www.gnu.org/software/dap/dap.html A FAIRE : DAP se voulait initialement un équivalent libre de SAS, mais il en est vraiment très loin. C'est plutôt une bibliothèque pour faire des calculs statistiques en C. A FAIRE : mettre une liste de bibliothèques (en C, C++, voire en Java) pour faire des calculs statistiques.
Un autre problème de R, c'est qu'il charge toutes les données qu'il manipule en mémoire : pour des données très volumineuses, ce n'est pas une bonne idée. (Toutefois, on peut accéder à un SGBD depuis R, et certains (PostgreSQL) acceptent des procédures stockées en R.) J'ai rencontré le problème une seule fois, en voulant manipuler un gros morceau de musique (quelques minutes) comme une série temporelle.
Il permet de créer des graphiques honnètes (sur ce point, il surpasse S-Plus), mais n'est pas pour autant un logiciel de présentation de données : pour des dessins informatifs, c'est bien, pour des brochures publicitaires, c'est insuffisant -- mais c'est une bonne base : on peut partir d'un dessin fait avec R pour l'améliorer avec d'autres logiciels -- ainsi, l'image de titre de ce document a réalisée avec R, PoVRay et Gimp.
Enfin, il s'agit d'un logiciel libre (attention, certaines de ses bibliothèques ne sont pas libres, mais simplement utilisables librement à des fins non commerciales -- mais la plupart sont libres).
A FAIRE : donner quelques exemples, si possible impressionnants.
Le document suivant explique comment utiliser R, progressivement, à un niveau élémentaire (il ne faut pas être effrayé par les expression « écart-type » ou « distribution gaussienne », mais c'est tout).
http://cran.r-project.org/doc/contrib/SimpleR.pdf
Un autre document élémentaire, avec des exercices, pour les germanophones :
http://cran.r-project.org/doc/contrib/s.pdf
Le manuel de référence de R est disponible directement sous R : il suffit de taper le nom de la commande sur laquelle on veut des éclaircissements précédé d'un point d'interrogation.
?round
Pour les mots-clefs du langage, ou pour des commandes non alphanumériques, il faut mettre des guillemets.
?"for" ?"[["
Si on ne connait pas le nom de la commande :
apropos("stem")
help.search("stem")La documentation en HTML (R lance konqueror) :
help.start()
Il y a aussi quelques fichiers PDF dans /usr/lib/R/library/*/doc/ : il s'agit parfois simplement d'une version PDF du manuel (quand le nom du fichier est celui de la bibliothèque), parfois d'explications différentes, plus pédagogiques.
% cd /usr/lib/R/library/ % ls */doc/*pdf CGIwithR/doc/CGIwithR-overview.pdf curve/doc/curve-guide.pdf DBI/doc/DBI.pdf dr/doc/drdoc.pdf dse1/doc/dse1-guide-012.pdf dse1/doc/dse1-guide.pdf dse1/doc/dse-guide.pdf dse2/doc/dse2-guide-012.pdf dse2/doc/dse2-guide-022.pdf dse2/doc/dse2-guide-024.pdf dse2/doc/dse2-guide-034.pdf dse2/doc/dse2-guide-036.pdf dse2/doc/dse2-guide-038.pdf dse2/doc/dse2-guide-040.pdf dse2/doc/dse2-guide-044.pdf dse2/doc/dse2-guide.pdf dsepadi/doc/dsepadi-guide.pdf e1071/doc/svmdoc.pdf evd/doc/guide.pdf genetics/doc/genetics_article.pdf geoRglm/doc/bookchap.pdf geoRglm/doc/geoRglm.intro.pdf gregmisc/doc/gregmisc.pdf grid/doc/uguide_0.7-4.pdf hwde/doc/hwde.pdf ifs/doc/IacLat2001.pdf ifs/doc/IacLat2002.pdf ipred/doc/ipred-examples.pdf juice/doc/juice-guide.pdf lasso2/doc/Manual-wo-help.pdf lmtest/doc/lmtest-intro.pdf maxstat/doc/maxstat.pdf monitor/doc/monitor-guide.pdf msm/doc/msm-manual.pdf multcomp/doc/Rmc.pdf permax/doc/Permax_Computations.pdf RSQLite/doc/DBI.pdf spatstat/doc/Intro.pdf spatstat/doc/Quickref.pdf StatDataML/doc/StatDataML.pdf strucchange/doc/strucchange-intro.pdf syskern/doc/syskern-guide.pdf tframe/doc/tframe-guide.pdf
Le document suivant (très instructif, même s'il est écrit dans une langue que je ne connais pas -- probablement de l'espagnol), explique en détails quels types de graphiques on peut réaliser avec R. J'en conseille vivement la lecture.
http://cran.r-project.org/doc/contrib/grafi3.pdf
L'anova en psychologie :
http://cran.r-project.org/doc/contrib/rpsych.html
La régression linéaire :
http://cran.r-project.org/doc/contrib/Faraway-PRA.pdf
Tout (c'est relativement complet, mais très dense : sans une connaissance préalable des sujets, c'est incompréhensible, mais sinon c'est très bien -- il y a même des exercices ) :
http://cran.r-project.org/doc/contrib/usingR.pdf
Le manuel est disponible en HTML (cherchez sur le disque, dans /usr/lib/R/), mais il n'y a pas les illustrations... On peut les rajouter comme suit (c'est un peu violent).
Pour que ce soit complet et à jour, j'installe la dernière version de R.
urpmi readline-devel wget http://stat.ethz.ch/CRAN/src/base/R-1.6.2.tgz tar zxvf R-1.6.2.tgz cd R-1.6.2/ ./configure --prefix=$HOME/spool/JUNK/R --with-readline --enable-R-shlib make -j 5 make install PATH=$HOME/spool/JUNK/R/bin:$PATH
avec tous les paquetages.
cd .. wget -r -l 0 -np -nc -x -k -E -U Mozilla -p -Q 50M cran.r-project.org/src/contrib/ for i in cran.r-project.org/src/contrib/*.tar.gz do echo Installing $i >> nohup.out echo Installing $i nohup R CMD INSTALL $i done
Ca plante pour les paquetages suivants : hdf5, Matrix, mimR, netCDF, Rmpi, RmSQL, RMySQL, RODBC, ROracle, rpvm, RQuantLib, rsprng, RSQLite, subselect, XML -- soit 15 erreurs pour 173 réussites. Très souvent, le problème vient de l'absence de certains paquetages de développement (comme readline, ci-dessus) : toutes les bibliothèques requises sont là, mais il y a juste les fichiers nécessaires à l'utilisation de programmes déjà compilés les requérant (les fichiers *.so) et il manque les fichiers *.h.
On peut ensuite récupérer tous les exemples,
#! perl -w
use strict;
my $n = 0;
# Ecriture des fichiers R
mkdir "Rdoc" || die "Cannot mkdir Rdoc/: $!";
my @libraries= `ls lib/R/library/`;
foreach my $lib (@libraries) {
chomp($lib);
print STDERR "Processing library \"$lib\"\n";
print STDERR `pwd`;
my @pages = grep /\.R$/, `ls lib/R/library/$lib/R-ex/`;
chdir "Rdoc" || die "Cannot chdir to Rdoc: $!";
mkdir "$lib" || die "Cannot mkdir $lib: $!";
chdir "$lib" || die "Cannot chdir to $lib: $!";
open(M, '>', "Makefile") || die "Cannot open Makefile for writing: $!";
print M "all:\n";
foreach my $page (@pages) {
chomp($page);
print STDERR " Processing man page \"$page\" in library \"$lib\"\n";
my $res = "";
$res .= "library($lib)\n";
$res .= "library(ts)##VZ##\n";
$res .= "library(lattice)##VZ##\n";
$res .= "library(nls)##VZ##\n";
$res .= "library(nlme)##VZ##\n";
$res .= "library(MASS)##VZ##\n";
$res .= "identify <- function (...) {}##VZ##\n";
$res .= "x11()##VZ##\n";
open(P, '<', "../../lib/R/library/$lib/R-ex/$page") ||
die "Cannot open lib/R/library/$lib/R-ex/$page for reading: $!";
# On marque les lignes où il faut copier l'écran (entre deux commandes)
while(<P>) {
s/^([^ #}])/try(invisible(dev.print(png,width=600,height=600,bg="white",filename="doc$n.png")))##VZ##\n$1/
&& $n++;
$res .= $_;
}
$res .= "try(invisible(dev.print(png,width=600,height=600,bg=\"white\",filename=\"doc$n.png\")))##VZ##\n"; $n++;
close P;
# On enlève la ligne supplémentaire dans les cas suivants :
# La ligne précdente se termine par une virgule, une parenthèse ouvrante, un signe "égale"
$res =~ s/[,(=+]\s*\n.*##VZ##.*?\n/,\n/g;
# La ligne précédente est vide
$res =~ s/^\s*\n.*##VZ##.*\n/\n/gm;
# La ligne précédente ne contient qu'un commentaire
$res =~ s/^(\s*#.*)\n.*##VZ##.*\n/$1\n/gm;
# La ligne suivante commence par une accolade ouvrante A FAIRE : vérifier (boot / abc.ci)
$res =~ s/^.*##VZ##.*\n\s*\{/\{/gm;
# On écrit le numéro correctement
$res =~ s/doc([0-9]).png/doc00000$1.png/g;
$res =~ s/doc([0-9][0-9]).png/doc0000$1.png/g;
$res =~ s/doc([0-9][0-9][0-9]).png/doc000$1.png/g;
$res =~ s/doc([0-9][0-9][0-9][0-9]).png/doc00$1.png/g;
$res =~ s/doc([0-9][0-9][0-9][0-9][0-9]).png/doc0$1.png/g;
open(W, ">", "${lib}_${page}") || die "Cannot open ${lib}_${page} for writing: $!";
print W $res;
close W;
print M "\tR --vanilla <${lib}_${page} >${lib}_${page}.out\n";
my $p = $page;
$p =~ s/\.R$//;
system 'cp', "../../lib/R/library/$lib/html/$p.html", "${lib}_$p.html";
}
print M "\ttouch all\n";
close(M);
chdir "../../" || die "Cannot chdir to ../../: $!";
}on les compile (compter quelques heures et quelques plantages),
cd Rdoc/
for i in *
do
(
cd $i
make
)
doneJe préfère donc paralléliser tout ça :
killall xscreensaver \ls -d */ | perl -p -e 's/(.*)/cd $1; make/' | perl fork.pl 5
où fork.pl permet de lancer plusieurs processus en même temps, mais pas trop :
#! /share/nfs/users1/umr-tge/zoonek/gnu/Linux/bin/perl -w
use strict;
my $MAX_PROCESSES = shift || 10;
use Parallel::ForkManager;
my $pm = new Parallel::ForkManager($MAX_PROCESSES);
while(<>){
my $pid = $pm->start and next;
system($_);
$pm->finish; # Terminates the child process
}On fait ensuite le ménage dans les fichiers PNG ainsi créés, et on écrit les fichiers HTML,
for i in */
do
(
cd $i
perl ../do.it_2.pl
)
doneoù le script do.it_2.pl est :
#! perl -w
use strict;
# Effaçons les fichiers PNG vides ou en double
print STDERR "Computing checksums\n";
use Digest::MD5 qw(md5);
my %checksum;
foreach my $f (sort(<*.png>)) {
if( -z $f ) {
unlink $f;
next;
}
local $/;
open(F, '<', $f) || warn "Cannot open $f for reading: $!";
my $m = md5(<F>);
close F;
if( exists $checksum{$m} ){
unlink $f;
} else {
$checksum{$m} = $f;
}
}
# Transformons tout cela en HTML...
print STDERR "Converting to HTML\n";
open(HTML, '>', "ALL.html") || warn "Cannot open ALL.html for writing: *!";
select(HTML);
print "<html>\n";
print "<head><title>R</title></head>\n";
print "<body>\n";
foreach my $f (<*.R.out>) {
my $page = $f;
$page =~ s/\.R.out$//;
# Lecture du fichier HTML initial
if( open(H, '<', "$page.html") ){
my $a = join '', <H>;
close H;
$a =~ s#^.*?<body>##gsi;
$a =~ s#<hr>.*?$##gsi;
print $a;
} else {
warn "Cannot open $page.html for reading: $!";
}
open(F, '<', $f) || warn "Cannot open $f for reading: $!";
#print "<h1>$f</h1>\n";
print "<h2>Worked out examples</h2>\n";
print "<pre>\n";
my $header=1;
while(<F>) {
if($header) {
$header=0 if m/to quit R/;
next;
}
if( m/(doc.*png).*##VZ##/ ){
my $png = $1;
next unless -f $png;
print "</pre>\n";
print "<img width=600 height=600 src=\"$png\">\n";
print "<pre>\n";
}
next if m/##VZ##/;
next if m/^>\s*###---/;
next if m/^>\s*##\s*___/;
next if m/^>\s*##\s*(alias|keywords)/i;
s/\&/\&/g;
s/</</g;
print;
}
close F;
print "</pre>\n";
print "<hr>\n";
}
print "</body>\n";
print "</html>\n";
close HTML;Pour une raison que j'ignore, les fichiers PNG ont un fond transparent : je le transforme en fond blanc avec ImageMagick (le fichier white.png est un fichier PNG tout blanc, de la même taille, 600x600, créé avec Gimp) -- là aussi, c'est très long...
for i in */*png do echo $i composite $i white.png $i done
Je n'ai pas tenu compte des éventuels liens dans les fichiers HTML : je les efface (ce n'est pas propre, mais...).
perl -p -i.bak -e 's#<a\s.*?>##gi; s#</a>##gi' **/ALL.html
Les fichiers HTML devraient plutôt s'appeler index.html :
rename 's/ALL.html/index.html/' */ALL.html
où rename est le programme :
#!/usr/bin/perl -w
use strict;
my $reg = shift @ARGV;
foreach (@ARGV) {
my $old = $_;
eval "$reg;";
die $@ if $@;
rename $old, $_;
}Voici un tout petit morceau du résultat (2.7Mo pour 268 images, alors que la documentation complète fait 93Mo pour 2078 images -- j'ai retiré une grande partie du code HTML)
Rdoc/index.html
En particulier, on pourra regarder les paquetages dont la documentation comporte le plus de graphiques.
% find -name "*.png" | perl -p -e 's#^\./##;s#/.*##' | sort | uniq -c | sort -n | tail -20 26 sm 27 cluster 29 ade4 29 nls 29 spatial 31 mgcv 34 car 38 cobs 41 MPV 44 mclust 45 MASS 45 vcd 49 grid 60 gregmisc 64 pastecs 70 qtl 77 splancs 88 strucchange 90 waveslim 113 spdep
C'est le mode statistique d'Emacs (sur la Mandrake, il n'est pas présent si on utilise Emacs, mais il l'est si on utilise XEmacs). On peut alors éditer du code en bénéficiant de l'indentation et de la colorisation automatique (emacs reconnait ces fichiers à leur extension),
on peut même lancer R sous Emacs (M-x R).
Dans cette partie, nous présentons les fonctions les plus simples de R, permettant de lire des données (ou de les construire, par simulation) et de les explorer numériquement ou graphiquement.
Paramètres statistiques et calculs :
> # 20 nombres entre 0 et 20,
> # arondis à un chiffre après la virgule
> x <- round(runif(20,0,20), digits=1)
> x
[1] 10.0 1.6 2.5 15.2 3.1 12.6 19.4 6.1 9.2 10.9 9.5 14.1 14.3 14.3 12.8
[16] 15.9 0.1 13.1 8.5 8.7
> min(x)
[1] 0.1
> max(x)
[1] 19.4
> median(x) # médiane
[1] 10.45
> mean(x) # moyenne
[1] 10.095
> var(x) # variance
[1] 27.43734
> sd(x) # standard deviation (écrat-type)
[1] 5.238067
> sqrt(var(x))
[1] 5.238067
> rank(x) # classement
[1] 10.0 2.0 3.0 18.0 4.0 12.0 20.0 5.0 8.0 11.0 9.0 15.0 16.5 16.5 13.0
[16] 19.0 1.0 14.0 6.0 7.0
> sum(x)
[1] 201.9
> length(x)
[1] 20
> round(x)
[1] 10 2 2 15 3 13 19 6 9 11 10 14 14 14 13 16 0 13 8 9
> fivenum(x) # les quantiles
[1] 0.10 7.30 10.45 14.20 19.40
> quantile(x) # les quantiles, calculés avec une convention différente
0% 25% 50% 75% 100%
0.10 7.90 10.45 14.15 19.40
> quantile(x, c(0,.33,.66,1))
0% 33% 66% 100%
0.100 8.835 12.962 19.400
> mad(x) # déviation moyenne à la médiane (normalisée)
[1] 5.55975
> cummax(x)
[1] 10.0 10.0 10.0 15.2 15.2 15.2 19.4 19.4 19.4 19.4 19.4 19.4 19.4 19.4 19.4
[16] 19.4 19.4 19.4 19.4 19.4
> cummin(x)
[1] 10.0 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6
[16] 1.6 0.1 0.1 0.1 0.1
> cor(x,sin(x/20)) # corrélation
[1] 0.997286Affichage d'un message ou d'une variable
print("Coucou !")Concaténation de chaines de caractères
paste("min: ", min (x$DS1, na.rm=T)))Ecriture dans un fichier (ou à l'écran)
cat("\\end{document}\n", file="RESULT.tex", append=TRUE)
Comme les logiciels de type matlab, R manipule des tableaux de nombres. Néanmoins, il y a différents types de tableaux : des vecteur (tableaux de dimension 1), des matrices (tableaux de dimension 2), des tableaux, des "Data Frames" (tableaux dans lesquels les colonnes ne sont pas nécessairement de même type), des listes (comme un tableau de dimension 1, mais les éléments ne sont pas nécessairement de même type). Nous allons maintenant présenter ces différents types de tableaux, les manières de les construire, des les manipuler, des les transformer, des les convertir.
Voici plusieurs manières de les définir.
> c(1,2,3,4,5) [1] 1 2 3 4 5 > 1:5 [1] 1 2 3 4 5 > seq(1,5,by=1) [1] 1 2 3 4 5 > seq(1,5,lenth=5) [1] 1 2 3 4 5
Voici quelques manières d'en récupérer un morceau.
> x <- seq(-1,1,by=.1) > x [1] -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 [16] 0.5 0.6 0.7 0.8 0.9 1.0 > x[5:10] [1] -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 > x[c(5,7:10)] [1] -0.6 -0.4 -0.3 -0.2 -0.1 > x[-(5:10)] # On enlève les éléments d'indice 5 à 10 [1] -1.0 -0.9 -0.8 -0.7 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 > x>0 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE [13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE > x[ x>0 ] [1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
On peut donner des noms aux composantes d'un vecteur.
> names(x) NULL > names(x) <- letters[1:length(x)] > x a b c d e f g h i j k l m n o p -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 q r s t u 0.6 0.7 0.8 0.9 1.0 > x["r"] r 0.7
On peut aussi définir ces noms lors de la création du vecteur.
> c(a=1, b=5, c=10, d=7) a b c d 1 5 10 7
Quelques opérations sur les vecteurs :
> x <- rnorm(10) > sort(x) [1] -1.4159893 -1.1159279 -1.0598020 -0.2314716 0.3117607 0.5376470 [7] 0.6922798 0.9316789 0.9761509 1.1022298 > rev(sort(x)) [1] 1.1022298 0.9761509 0.9316789 0.6922798 0.5376470 0.3117607 [7] -0.2314716 -1.0598020 -1.1159279 -1.4159893 > o <- order(x) > o [1] 3 1 9 6 4 7 8 10 2 5 > x[ o[1:3] ] [1] -1.415989 -1.115928 -1.059802 > x <- sample(1:5, 10, replace=T) > x [1] 1 4 5 3 1 3 4 5 3 1 > sort(x) [1] 1 1 1 3 3 3 4 4 5 5 > unique(x) # Il n'est pas nécassaire de trier les données [1] 1 4 5 3
Voici d'autres manières de créer des vecteurs. La commande seq permet de créer des vecteurs contenant les termes d'une suite arithmétique.
> seq(0,10, length=11) [1] 0 1 2 3 4 5 6 7 8 9 10 > seq(0,10, by=1) [1] 0 1 2 3 4 5 6 7 8 9 10
La commande rep permet de répéter un nombre ou un vecteur.
> rep(1,10) [1] 1 1 1 1 1 1 1 1 1 1 > rep(1:5,3) [1] 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5
On peut aussi lui demander de répéter chaque élément plusieurs fois.
> rep(1:5,each=3) [1] 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5
On peut mélanger les deux opérations précédentes.
> rep(1:5,2,each=3) [1] 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5
La commande gl fait des choses comparables, essentiellement pour créer des facteurs : nous en reparlerons plus loin.
Il s'agit de vecteurs codant une variable qualitative. On les crée à l'aide de la commande factor.
> x <- factor( sample(c("Oui", "Non", "Peut-être"), 5, replace=T) )
> x
[1] Peut-être Peut-être Peut-être Peut-être Non
Levels: Non Peut-êtreOn peut préciser la liste des valeurs possibles de cette variable qualitative.
> l <- c("Oui", "Non", "Peut-être")
> x <- factor( sample(l, 5, replace=T), levels=l )
> x
[1] Non Peut-être Non Oui Oui
Levels: Oui Non Peut-être
> levels(x)
[1] "Oui" "Non" "Peut-être"On peut préférer un tableau avec les effectifs.
> table(x)
x
Oui Non Peut-être
2 2 1Si les valeurs suivant un certain motif, on peut les créer à l'aide de la commande gl.
> gl(1,4) [1] 1 1 1 1 Levels: 1 > gl(2,4) [1] 1 1 1 1 2 2 2 2 Levels: 1 2 > gl(2,4, labels=c(T,F)) [1] TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE Levels: TRUE FALSE > gl(2,1,8) [1] 1 2 1 2 1 2 1 2 Levels: 1 2 > gl(2,1,8, labels=c(T,F)) [1] TRUE FALSE TRUE FALSE TRUE FALSE TRUE FALSE Levels: TRUE FALSE
On peut faire le produit cartésien de deux facteurs.
> x <- gl(2,4) > x [1] 1 1 1 1 2 2 2 2 Levels: 1 2 > y <- gl(2,1,8) > y [1] 1 2 1 2 1 2 1 2 Levels: 1 2 > interaction(x,y) [1] 1.1 1.2 1.1 1.2 2.1 2.2 2.1 2.2 Levels: 1.1 2.1 1.2 2.2
La commande expand.grid est comparable (mais elle produit un data.frame).
> x <- c("A", "B", "C")
> y <- 1:2
> z <- c("a", "b")
> expand.grid(x,y,z)
Var1 Var2 Var3
1 A 1 a
2 B 1 a
3 C 1 a
4 A 2 a
5 B 2 a
6 C 2 a
7 A 1 b
8 B 1 b
9 C 1 b
10 A 2 b
11 B 2 b
12 C 2 b
Il s'agit d'une "liste" de vecteurs, chacun ayant la même longueur. Généralement, les lignes du tableau sont les individus et les colonnes les différentes mesures sur ces individus.
> n <- 10 > df <- data.frame( x=rnorm(n), y=sample(c(T,F),n,replace=T) )
La commande str affiche la structure d'un objet (n'importe quel type d'objet) et une partie des données qu'il contient.
> str(df) `data.frame': 10 obs. of 2 variables: $ x: num 0.515 -1.174 -0.523 -0.146 0.410 ... $ y: logi FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
Quand on manipule des objets plus compliqués, qui ont leur propre commande print, la commande unclass peut s'avérer plus claire.
n <- 10 x <- runif(n) y <- 1 - 2 * x + rnorm(n) r <- lm(y~x) r # Appelle print.lm str(r) print.default(r) unclass(r) # Pareil print.lm # Pour avoir une idée de ce que l'on peut faire avec l'objet... summary.lm print.summary.lm
La commande summary affiche des informations assez compactes sur un objet (ici, un data.frame, mais ça marche avec presque tout).
> summary(df)
x y
Min. :-1.17351 Length:10
1st Qu.:-0.42901 Mode :logical
Median : 0.13737
Mean : 0.09217
3rd Qu.: 0.48867
Max. : 1.34213
> df
x y
1 0.51481130 FALSE
2 -1.17350867 TRUE
3 -0.52338041 FALSE
4 -0.14589347 FALSE
5 0.41022626 FALSE
6 1.34213009 TRUE
7 0.77715729 FALSE
8 -0.55460889 FALSE
9 -0.03843468 FALSE
10 0.31318467 FALSEOn peut accèder aux colonnes de plusieurs manières.
> df$x [1] 0.51481130 -1.17350867 -0.52338041 -0.14589347 0.41022626 1.34213009 [7] 0.77715729 -0.55460889 -0.03843468 0.31318467 > df[,1] [1] 0.51481130 -1.17350867 -0.52338041 -0.14589347 0.41022626 1.34213009 [7] 0.77715729 -0.55460889 -0.03843468 0.31318467 > df[["x"]] [1] 0.51481130 -1.17350867 -0.52338041 -0.14589347 0.41022626 1.34213009 [7] 0.77715729 -0.55460889 -0.03843468 0.31318467 > dim(df) [1] 10 2 > names(df) [1] "x" "y" > row.names(df) [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10"
On peut changer le nom des lignes ou des colonnes.
> names(df) <- c("a", "b")
> row.names(df) <- LETTERS[1:10]
> names(df)
[1] "a" "b"
> row.names(df)
[1] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J"
> str(df)
`data.frame': 10 obs. of 2 variables:
$ a: num 0.515 -1.174 -0.523 -0.146 0.410 ...
$ b: logi FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSEOn peut transformer les colonnes d'un Data Frame en variables à l'aide de la commade attach() (c'est le même principe que les espaces de nommage (namespaces) du C++). Il ne faut pas oublier de détacher le Data Frame quand on a fini de jouer avec, sinon les espaces de nommage risquent de s'empiler.
> data(faithful) > str(faithful) `data.frame': 272 obs. of 2 variables: $ eruptions: num 3.60 1.80 3.33 2.28 4.53 ... $ waiting : num 79 54 74 62 85 55 88 85 51 85 ... > attach(faithful) > str(eruptions) num [1:272] 3.60 1.80 3.33 2.28 4.53 ... > detach()
La commande merge permet de fusionner deux data.frames. C'est l'analogue du JOIN des bases de données. On dispose de deux data.frames a (de colonnes x, y, z) et b (de colonnes x1, x2, y, z) et certaines des observations (lignes) de a correspondent à certaines de celles de b : en les fusionnant, on peut obtenir un data.frame ayant pour colonnes x, x1, x2, y, z. Dans cet exemple, la commande
merge(a,b)
est équivelente à
SELECT * FROM a,b WHERE a.y = b.y AND a.z = b.z
Par défaut, la jonction se fait sur toutes les colonnes communes aux deux data.frames, mais on peut ne prendre que certaines d'entre elles, grace à l'argument by.
On utilise souvent les data.frames pour stocker des données à analyser. Nous détaillerons ces exemples un peu plus tard.
# Régression data(cars) lm( dist ~ speed, data=cars) # Régression polynomiale lm( dist ~ poly(speed,3), data=cars) # Régression logistique A FAIRE glm( ..., family=binomial, data=...) # Régression non linéaire A FAIRE # Analyse en composantes principales data(USArrest) princomp( ~ Murder + Assault + UrbanPop, data=USArrest) # Graphiques en treillis A FAIRE
Nous verrons un peu plus loin comment transformer les data.frames, car il y a plusieurs manières de mettre le résultat d'une expérience dans un tableau -- mais généralement, on préférera la présentation avec le plus de lignes et le moins de colonnes.
Dans des exemples simples et de petites tailles, on peut taper directement les données qu'on désire analyser. Dans ce document, nous utiliserons beaucoup de données simulées : elles sont plus volumineuses, mais il suffit de taper quelques lignes pour les obtenir.
A contrario, dans des situations réelles, les données sont volumineuses et stockées dans des fichiers ou des bases de données : comment les importer sous R ?
Personnellement, j'utilise souvent la commande source, qui n'est pourtant pas du tout faite pour ça : elle lit du code, pas des données. Les données que je devais analyser étaient sous une forme assez peu standard (des alignement multiples de séquences génomiques) : j'ai donc écrit un petit programme en Perl pour convertir ce format en code pour R.
Plus précisément, mes données ressemblaient à
CLUSTAL W (1.83) multiple sequence alignment
AB020905 ATGACCAACATCCGAAAAACCCACCCATTAGCTAAAATCATCAACAACTCATTTATTGAC
AB020906 ATGACCAACATCCGAAAAACCCACCCATTAGCTAAAATCATCAACAACTCATTTATTGAC
AB020907 ATGACCAACATCCGAAAAACCCACCCATTAGCTAAAATCATCAACAACTCACTTATTGAC
AB020908 ATGACCAACATCCGAAAAACCCACCCATTAGCTAAAATCATCAACAACTCATTTATTGAC
AB020909 ATGACCAACATCCGAAAAACCCACCCATTAGCTAAAATCATCAACAACTCATTTATTGAC
*************************************************** ********
AB020905 CTTCCAACACCATCAAACATCTCGGCATGATGAAACTTTGGATCCCTCCTTGGAGTATGT
AB020906 CTTCCAACACCATCAAACATCTCAGCATGATGAAACTTTGGATCCCTCCTTGGAGTATGT
AB020907 CTTCCAACACCATCAAACATCTCAGCATGATGAAACTTTGGATCCCTCCTTGGAGTATGT
AB020908 CTTCCAACACCATCAAACATCTCAGCATGATGAAACTTTGGATCCCTCCTCGGAGTATGT
AB020909 CTTCCAACACCATCAAACATCTCAGCATGATGAAACTTTGGATCCCTCCTCGGAGTATGT
*********************** ************************** *********le programme était
#! perl -w
use strict;
my @seq;
my @names;
my $i=0;
# On se place après la première ligne vide
while (<>) {
chomp;
print STDERR "Skipping $. ($_)\n";
last if m/^\s*$/;
}
while (<>) {
chomp;
if( m/^\s*$/ ){
$i=0;
print STDERR "Skipping $. ($_)\n";
next;
}
print STDERR "Reading $. ($i) ($_)\n";
if (m/^([^\s]+?)\s+(.*)/) {
print STDERR "Remembering $.\n";
$names[$i] = $1;
$seq[$i] .= $2;
}
$i++;
}
# foreach my $s (@seq) { print "$s\n"; }
print "d <- matrix( c(\n";
foreach my $s (@seq) {
print '"'. join('", "', split('', $s)) .'",' ."\n";
}
print "), nr=". (scalar @seq) .", byrow=T)\n";
print "rownames(d) <- c('". join("', '", @names) ."')\n"et le résultat ressemblait à
d <- matrix( c( "A", "T", "G", "A", "C", "C", "A", "A", "C", "A",
"T", "C", "C", "G", "A", "A", "A", "A", "A", "C", "C", "C", "A",
...
), nr=5, byrow=T)
rownames(d) <- c('AB020905', 'AB020906', 'AB020907', 'AB020908', 'AB020909')Le problème, avec cette méthode, c'est que ce ne sont pas des données, mais du code. Si on veut juste les utiliser avec R et pas avec un autre logiciel, c'est bon, mais sinon, mieux vaut un format un peu plus universel.
La commande read.table permet de lire des data.frames, i.e., des tableaux (rectangulaires), dont les colonnes peuvent éventuellement avoir des types différents (mais le type ne varie pas au sein d'une colonne). Si on reprend l'exemple précédent, il s'agirait simplement de produire un fichier qui ressemble à
AB020905 "T" "T" "A" "A" "A" "G" "T" "G" ... AB020906 "T" "T" "A" "A" "A" "G" "T" "G" ... AB020907 "T" "T" "A" "A" "A" "G" "T" "G" ... AB020908 "T" "T" "A" "A" "A" "G" "T" "G" ... AB020909 "T" "T" "A" "A" "T" "G" "T" "G" ...
(avec des lignes très, très longues).
Souvent, ça marche sans qu'on ait à y penser, mais des problèmes peuvent survenir.
Considérons tout d'abord le cas simple où le fichier ne contient que les données, des données numériques, sans nom de ligne ou de colonne. Il pourrait ressembler à
2 7 3 9 2 8 7 3 2 2 6 2 8 8 1
Il vient :
> d <- read.table('A.txt')
> d
V1 V2 V3 V4 V5
1 2 7 3 9 2
2 8 7 3 2 2
3 6 2 8 8 1On remarquera que R a donné des noms aux colonnes. S'ils ne nous plaisent pas, on peut les changer.
> names(d) [1] "V1" "V2" "V3" "V4" "V5" > length(d) [1] 5 > names(d) <- 1:length(d) > d 1 2 3 4 5 1 2 7 3 9 2 2 8 7 3 2 2 3 6 2 8 8 1 > names(d) <- LETTERS[1:length(d)] > d A B C D E 1 2 7 3 9 2 2 8 7 3 2 2 3 6 2 8 8 1
Notre fichier peut être plus complexe et contenir des nom de lignes.
x1 2 7 3 9 2 x2 8 7 3 2 2 x3 6 2 8 8 1
Mais ça marche un peu moins bien, car l'ordinateur n'a aucun moyen de savoir que la première colonne est le nom des lignes et pas une variable.
> read.table('A.txt')
V1 V2 V3 V4 V5 V6
1 x1 2 7 3 9 2
2 x2 8 7 3 2 2
3 x3 6 2 8 8 1On peut lui demander de retirer la première colonne et de l'utiliser comme nom de lignes (c'est un bon exercice : essayez de le faire seul)
> d <- read.table('A.txt')
> row.names(d) <- d[,1]
> d <- d[,-1]
> d
V2 V3 V4 V5 V6
x1 2 7 3 9 2
x2 8 7 3 2 2
x3 6 2 8 8 1
> names(d) <- LETTERS[1:length(d)]
> d
A B C D E
x1 2 7 3 9 2
x2 8 7 3 2 2
x3 6 2 8 8 1Autre situation : on a à la fois des noms de lignes et des noms de colonnes. Le fichier ressemble à :
A B C D E x1 2 7 3 9 2 x2 8 7 3 2 2 x3 6 2 8 8 1
Maintenant, R s'apperçoit que la première ligne contient le nom des variables et que la première colonne contient le nom des observations, car la première ligne du fichier comporte un champ de moins que les autres lignes.
> read.table('A')
A B C D E
x1 2 7 3 9 2
x2 8 7 3 2 2
x3 6 2 8 8 1Dernière situation : les colonnes ont des noms, mais pas les lignes. Le fichier ressemble à
A B C D E 2 7 3 9 2 8 7 3 2 2 6 2 8 8 1
Si on essaye naïvement :
> d <- read.table('A.txt')
> d
V1 V2 V3 V4 V5
1 A B C D E
2 2 7 3 9 2
3 8 7 3 2 2
4 6 2 8 8 1
> str(d)
`data.frame': 4 obs. of 5 variables:
$ V1: Factor w/ 4 levels "2","6","8","A": 4 1 3 2
$ V2: Factor w/ 3 levels "2","7","B": 3 2 2 1
$ V3: Factor w/ 3 levels "3","8","C": 3 1 1 2
$ V4: Factor w/ 4 levels "2","8","9","D": 4 3 1 2
$ V5: Factor w/ 3 levels "1","2","E": 3 2 2 1D'une part, il ne s'est pas apperçu que la première ligne contenait le nom des colonnes, d'autre part, il a cru que chaque colonne contenait des chaines de caractères... on peut contourner ces problèmes en ajoutant un argument à la commande read.table.
> read.table('A.txt', header=T)
A B C D E
1 2 7 3 9 2
2 8 7 3 2 2
3 6 2 8 8 1Avons-nous fait le tour de l'utilisation et des problèmes soulevés par la commande read.table ? Eh bien non. Reprenons le tout premier exemple, avec les séquences nucléiques.
> read.table('A.txt')
V1 V2 V3 V4 V5 V6 V7 V8 V9
1 AB020905 TRUE TRUE A A A G TRUE G
2 AB020906 TRUE TRUE A A A G TRUE G
3 AB020907 TRUE TRUE A A A G TRUE G
4 AB020908 TRUE TRUE A A A G TRUE G
5 AB020909 TRUE TRUE A A T G TRUE GChaque colonne contenait des caractères (les quatres lettres A, C, G, T), mais il a compris que les T étaient des valeurs booléennes. On peut heureusement préciser le type des colonnes (ici, c'est toujours le même type, donc il suffit de le mentionner une seule fois).
> read.table('A', colClasses=c('character'))
V1 V2 V3 V4 V5 V6 V7 V8 V9
1 AB020905 T T A A A G T G
2 AB020906 T T A A A G T G
3 AB020907 T T A A A G T G
4 AB020908 T T A A A G T G
5 AB020909 T T A A T G T GA FAIRE : et si le fichier n'est pas sur le disque local, mais donné par un URL ?
A FAIRE : et si le fichier provient d'un tableur (dont je taierai le nom...)
A FAIRE : la commande scan
A FAIRE : connection à des bases de données.
A FAIRE : lire la documentation sur l'importation de données sous R.
Les vecteurs contiennent uniquement des types simples (nombres, booléens ou chaînes de caractères) ; les listes, au contraire, peuvent contenir n'importe quoi, par exemples des Data Frames ou d'autres listes. On peut les utiliser pour stocker des données complexes, comme des arbres. On peut aussi les utiliser, simplement, comme des tables de hachage.
> h <- list()
> h[["foo"]] <- 1
> h[["bar"]] <- c("a", "b", "c")
> str(h)
List of 2
$ foo: num 1
$ bar: chr [1:3] "a" "b" "c"Par exemple, les paramètres graphiques sont stockés dans une liste utilisée comme table de hachage.
> str( par() ) List of 68 $ adj : num 0.5 $ ann : logi TRUE $ ask : logi FALSE $ bg : chr "transparent" $ bty : chr "o" $ cex : num 1 $ cex.axis : num 1 $ cex.lab : num 1 $ cex.main : num 1.2 $ cex.sub : num 1 $ cin : num [1:2] 0.147 0.200 ... $ xpd : logi FALSE $ yaxp : num [1:3] 0 1 5 $ yaxs : chr "r" $ yaxt : chr "s" $ ylog : logi FALSE
Pour effacer un élément d'une liste :
> h[["bar"]] <- NULL > str(h) List of 1 $ foo: num 1
Les matrices sont des tableaux, mais contrairement aux data.frames (dont le type peut carier d'une colonne à l'autre), leurs éléments sont tous du même type.
Une matrice :
> m <- matrix( c(1,2,3,4), nrow=2 )
> m
[,1] [,2]
[1,] 1 3
[2,] 2 4ATTENTION : par défaut, les coefficients des matrices sont écrits verticalement.
> matrix( 1:3, nrow=3, ncol=3 )
[,1] [,2] [,3]
[1,] 1 1 1
[2,] 2 2 2
[3,] 3 3 3
> matrix( 1:3, nrow=3, ncol=3, byrow=T )
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 1 2 3
[3,] 1 2 3
> t(matrix( 1:3, nrow=3, ncol=3 ))
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 1 2 3
[3,] 1 2 3Le produit matriciel :
> x <- matrix( c(6,7), nrow=2 )
> x
[,1]
[1,] 6
[2,] 7
> m %*% x
[,1]
[1,] 27
[2,] 40Le déterminant d'une matrice :
> det(m) [1] -2
La transposée d'une matrice :
> t(m)
[,1] [,2]
[1,] 1 2
[2,] 3 4Une matrice diagonale :
> diag(c(1,2))
[,1] [,2]
[1,] 1 0
[2,] 0 2La matrice identité (ou la matrice d'une homothétie) :
> diag(1,2)
[,1] [,2]
[1,] 1 0
[2,] 0 1
> diag(rep(1,2))
[,1] [,2]
[1,] 1 0
[2,] 0 1
> diag(2)
[,1] [,2]
[1,] 1 0
[2,] 0 1On a déjà vu les fonctions cbind et rbind, qui permettent de juxtaposer des vecteurs-colonnes ou d'empiler des vecteurs-lignes pour créer des matrices. On peut aussi les utiliser avec des matrices.
> cbind( c(1,2), c(3,4) )
[,1] [,2]
[1,] 1 3
[2,] 2 4
> rbind( c(1,3), c(2,4) )
[,1] [,2]
[1,] 1 3
[2,] 2 4La trace d'une matrice :
> sum(diag(m)) [1] 5
L'inverse d'une matrice :
> solve(m)
[,1] [,2]
[1,] -2 1.5
[2,] 1 -0.5Généralement, on n'a pas besoin de l'inverse d'une matrice, mais simplement de multiplier par l'inverse d'une matrice : c'est plus rapide et numériquement plus stable.
> solve(m, x)
[,1]
[1,] -1.5
[2,] 2.5
> solve(m) %*% x
[,1]
[1,] -1.5
[2,] 2.5Valeurs propres :
> eigen(m)$values [1] 5.3722813 -0.3722813
Vecteurs propres :
> eigen(m)$vectors
[,1] [,2]
[1,] -0.5742757 -0.9093767
[2,] -0.8369650 0.4159736On vérifie qu'on a effectivement diagonalisé la matrice :
> p <- eigen(m)$vectors
> d <- diag(eigen(m)$values)
> p %*% d %*% solve(p)
[,1] [,2]
[1,] 1 3
[2,] 2 4C'est peut-être le bon moment pour rappeler les principales décomposition d'une matrice, que l'on rencontre en algèbre linéaire.
La décomposition LU, ou plus précisément PA = LDU (P: matrice de permutation ; L, U : matrices triangulaires inférieure ou supérieure, avec des 1 sur la diagonale ; L contient toutes les opérations effectuées sur les lignes ; D : matrice des pivots) traduit le résultat de l'algorithme du pivot de Gauss.
On n'en a pas réellement besoin, car l'algorithme du pivot de Gauss est implémenté dans la commande solve, qui permet de multiplier par l'inverse d'une matrice (sans calculer explicitement cet inverse, ce qui serait numériquement plus instable).
?solve
La décomposition de Cholesky est un cas particulier de la décomposition LU : si A est symétrique réelle définie positive, alors on peut l'écrire A*A' où A est triangulaire supérieure. On l'utilise pour résoudre des systèmes d'équations AX=Y dans le cas où A est symétrique -- c'est le cas des équations définissant les estimateurs des moindres carrés. Dans le même ordre d'idées, on peut l'utiliser pour
Nous verrons une autre application pour simuler des variables normales non indépendantes, de matrice de covariance donnée.
La décomposition de Jordan est une généralisation de la diagonalisation (car toutes les matrices ne sont pas diagonalisables).
On n'en a pas réellement besoin, car si on prend une matrice au hasard, elle est diagonalisable.
Il y a ensuite plein de décompositions basées sur la matrice t(A)*A.
La décomposition A=QR (R : triangulaire supérieure, Q : unitaire) traduit l'orthonormalisation de Gram-Schmidt des colonnes de A (on peut l'obtenir à l'aide de la décomposition LU de t(A)*A).
?qr
C'est elle qu'on utilise généralement pour faire une régression :
Modèle : Y = b X + bruit X = QR \hat Y = Q Q' Y b = R^-1 Q' Y
La décomposition en valeurs singulières A=Q1*S*Q2 (Q1, Q2 : matrices des vecteurs propres de A*t(A) et t(A)*A ; S : matrice diagonale des racines carrées des valeurs propres de A*t(A) ou t(A)*A), qui donne, dans le cas où A est symétrique, sa diagonalisation dans une base orthonormale, et qui intervient aussi dans le calcul du pseudo-inverse. On peut aussi voir cette décomposition comme une écriture de la matrice comme somme de matrices de rang 1, de sorte que les premières approchent « au mieux » la matrice de départ.
?svd
La décomposition polaire, A = QR (Q : orthogonale, R : symétrique définie positive), qui est un analogue de la décomposition polaire d'un nombre complexe, et qui sépare la rotation des étirements. On la rencontre dans la méthode des moindres carrés : on cherche à minimiser la valeur absolue de Ax - b, on montre que cela se ramène à résoudre
t(A) A x = t(A) b
(généralement, t(A)A est inversible, sinon, on utiliserait des pseudo-inverses).
Il y a aussi un type "array", qui généralise les matrices en dimension quelconque.
> d <- array(rnorm(3*3*2), dim=c(3,3,2))
> d
, , 1
[,1] [,2] [,3]
[1,] 0.97323599 -0.7319138 -0.7355852
[2,] 0.06624588 -0.5732781 -0.4133584
[3,] 1.65808464 -1.3011671 -0.4556735
, , 2
[,1] [,2] [,3]
[1,] 0.6314685 0.6263645 1.2429024
[2,] -0.2562622 -1.5338054 0.9634999
[3,] 0.1652014 -0.9791350 -0.2040375
> str(d)
num [1:3, 1:3, 1:2] 0.9732 0.0662 1.6581 -0.7319 -0.5733 ...On les rencontre, par exemple, pour décrire des tableaux de contingence avec plus de deux variables.
> data(HairEyeColor)
> HairEyeColor
, , Sex = Male
Eye
Hair Brown Blue Hazel Green
Black 32 11 10 3
Brown 38 50 25 15
Red 10 10 7 7
Blond 3 30 5 8
, , Sex = Female
Eye
Hair Brown Blue Hazel Green
Black 36 9 5 2
Brown 81 34 29 14
Red 16 7 7 7
Blond 4 64 5 8
> str(HairEyeColor)
table [1:4, 1:4, 1:2] 32 38 10 3 11 50 10 30 10 25 ...
- attr(*, "dimnames")=List of 3
..$ Hair: chr [1:4] "Black" "Brown" "Red" "Blond"
..$ Eye : chr [1:4] "Brown" "Blue" "Hazel" "Green"
..$ Sex : chr [1:2] "Male" "Female"
- attr(*, "class")= chr "table"Le fait que la commande "str" nous dise "table" et pas "array" n'est pas inquétant :
> is.array(HairEyeColor) [1] TRUE
A FAIRE On peut attacher des méta-informations à un objet. Exemples
Nous allons voir tout de suite une autre application de l'un de ces attributs : la notion de classe.
L'affichage du résultat de certaines commandes ne ressemble pas du tout à celui des types que nous venons de décrire. C'est le cas par exemple de la régression
n <- 200 x <- rnorm(n) y <- 1 - 2 * x + rnorm(n) r1 <- lm(y~x) r2 <- summary(r1)
qui donne
> r1
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept) x
0.924 -2.042
> r2
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-2.85364 -0.66754 -0.04169 0.61238 2.78004
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.92395 0.07345 12.58 <2e-16 ***
x -2.04152 0.07613 -26.82 <2e-16 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 1.038 on 198 degrees of freedom
Multiple R-Squared: 0.7841, Adjusted R-squared: 0.783
F-statistic: 719.1 on 1 and 198 DF, p-value: < 2.2e-16Pourtant, la commande str nous dit qu'il s'agit quand-même de l'un des types que nous venons de décrire -- très souvent, comme ici, une liste (j'en ai enlevé un gros morceau).
> str(r2) List of 11 $ call : language lm(formula = y ~ x) $ terms :Classes 'terms', 'formula' length 3 y ~ x .. ..- attr(*, "variables")= language list(y, x) .. ..- attr(*, "factors")= int [1:2, 1] 0 1 .. .. ..- attr(*, "dimnames")=List of 2 .. .. .. ..$ : chr [1:2] "y" "x" .. .. .. ..$ : chr "x" (...) - attr(*, "class")= chr "summary.lm"
Il y a une différence notable avec les listes que nous avons créées précédemment : l'attribut "class". Les objets r1 et r2 sont ce classe "lm" et "summary.lm". cela a pour effet de changer certaines fonctions qui s'appliquent à ces objets : en particulier, la fonction print, qui est appelée implicitement quand on demande l'affichage de l'objet.
> print
function (x, ...)
UseMethod("print")
<environment: namespace:base>
> print.lm
function (x, digits = max(3, getOption("digits") - 3), ...)
{
cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
if (length(coef(x))) {
cat("Coefficients:\n")
print.default(format(coef(x), digits = digits), print.gap = 2,
quote = FALSE)
}
else cat("No coefficients\n")
cat("\n")
invisible(x)
}
<environment: namespace:base>On qualifie la commande print de "méthode" : son code sera différent selon la classe de l'objet auquel on l'applique.
On peut d'ailleurs écrire ses propres classes : il suffit d'ajouter un attribut à un objet et de définir la méthode correspondante.
x <- pi
attr(x, 'class') <- "nombre"
print.nombre <- function (x) {
cat("(nombre) ")
cat(signif(x))
cat("\n")
}Il vient :
> x (nombre) 3.14159
On peut aussi définir de nouvelles méthodes.
affiche <- function (x,...) {
UseMethod("affiche")
}
affiche.default <- print
affiche.nombre <- function (x) {
cat("(nombre) ")
cat(signif(x))
cat("\n")
}Il vient :
> affiche(x) (nombre) 3.14159 > affiche(pi) [1] 3.141593
La commande methods donne la liste des implémentations de cette méthode.
> methods(plot)
[1] plot.acf* plot.ACF* plot.augPred*
[4] plot.compareFits* plot.data.frame plot.decomposed.ts*
[7] plot.default plot.dendrogram* plot.density
[10] plot.factor plot.formula plot.function
[13] plot.gls* plot.hclust* plot.histogram
[16] plot.HoltWinters* plot.intervals.lmList* plot.isoreg*
[19] plot.lm plot.lme plot.lme1*
[22] plot.lmList* plot.mlm plot.nffGroupedData*
[25] plot.nfnGroupedData* plot.nls* plot.nmGroupedData*
[28] plot.pdMat* plot.POSIXct plot.POSIXlt
[31] plot.ppr* plot.prcomp* plot.princomp*
[34] plot.profile.nls* plot.ranef.lme* plot.ranef.lmList*
[37] plot.shingle* plot.simulate.lme* plot.spec
[40] plot.spec1* plot.spec.coherency plot.spec.phase
[43] plot.stl* plot.table plot.ts
[46] plot.tskernel* plot.TukeyHSD plot.Variogram*
Non-visible functions are asteriskedRemarquons qu'un objet peut avoir plusieurs types différents : l'attribut "class" peut contenur un vecteur de chaines de caractères -- lors de l'appel d'une méthode, on prend ces chaines une par une, jusqu'à trouver une méthode qui soit implémentée.
Par exemple, le résumtat de la commande aov
n <- 500 x <- rnorm(n) y <- 1 - x + rnorm(n) r <- aov(y~x)
a pour classe aov et lm.
> class(r) [1] "aov" "lm" > attr(r,"class") [1] "aov" "lm"
Si on veut utiliser un objet complexe, obtenu comme résultat d'une certaine commande, en en extrayant certains éléments, ou si on veut explorer son contenu, la fonction d'affichage ne suffit pas : il nous faut d'autres moyens de regarder à l'intérieur d'un objet.
La fonction unclass enlève la classe d'un objet : ainsi, son affichage est effectué par la méthode print.default et on voit réellement ce qu'il contient.
> data(USArrests)
> r <- princomp(USArrests)$loadings
> r
Loadings:
Comp.1 Comp.2 Comp.3 Comp.4
Murder 0.995
Assault -0.995
UrbanPop -0.977 -0.201
Rape -0.201 0.974
Comp.1 Comp.2 Comp.3 Comp.4
SS loadings 1.00 1.00 1.00 1.00
Proportion Var 0.25 0.25 0.25 0.25
Cumulative Var 0.25 0.50 0.75 1.00
> class(r)
[1] "loadings"
> unclass(r)
Comp.1 Comp.2 Comp.3 Comp.4
Murder -0.04170432 0.04482166 0.07989066 0.99492173
Assault -0.99522128 0.05876003 -0.06756974 -0.03893830
UrbanPop -0.04633575 -0.97685748 -0.20054629 0.05816914
Rape -0.07515550 -0.20071807 0.97408059 -0.07232502On peut aussi utiliser directement la fonction print.default().
> print.default(r)
Comp.1 Comp.2 Comp.3 Comp.4
Murder -0.04170432 0.04482166 0.07989066 0.99492173
Assault -0.99522128 0.05876003 -0.06756974 -0.03893830
UrbanPop -0.04633575 -0.97685748 -0.20054629 0.05816914
Rape -0.07515550 -0.20071807 0.97408059 -0.07232502
attr(,"class")
[1] "loadings"La fonction str affiche le contenu d'un objet en tronquant tous les vecteurs : cela permet de voir le contenu d'objets gros mais pas nécessairement complexes.
> str(r) loadings [1:4, 1:4] -0.0417 -0.9952 -0.0463 -0.0752 0.0448 ... - attr(*, "dimnames")=List of 2 ..$ : chr [1:4] "Murder" "Assault" "UrbanPop" "Rape" ..$ : chr [1:4] "Comp.1" "Comp.2" "Comp.3" "Comp.4" - attr(*, "class")= chr "loadings" > str(USArrests) `data.frame': 50 obs. of 4 variables: $ Murder : num 13.2 10 8.1 8.8 9 7.9 3.3 5.9 15.4 17.4 ... $ Assault : int 236 263 294 190 276 204 110 238 335 211 ... $ UrbanPop: int 58 48 80 50 91 78 77 72 80 60 ... $ Rape : num 21.2 44.5 31 19.5 40.6 38.7 11.1 15.8 31.9 25.8 ...
Enfin, si on veut avoir une idée de ce qu'on peut faire avec le contenu d'un objet, on peut toujours regarder le contenu des méthodes print ou summary qui s'y appliquent.
> print.lm
function (x, digits = max(3, getOption("digits") - 3), ...)
{
cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
if (length(coef(x))) {
cat("Coefficients:\n")
print.default(format(coef(x), digits = digits), print.gap = 2,
quote = FALSE)
}
else cat("No coefficients\n")
cat("\n")
invisible(x)
}
<environment: namespace:base>
> summary.lm
function (object, correlation = FALSE, symbolic.cor = FALSE,
...)
{
z <- object
p <- z$rank
if (p == 0) {
r <- z$residuals
n <- length(r)
etc.
> print.summary.lm
function (x, digits = max(3, getOption("digits") - 3), symbolic.cor = x$symbolic.cor,
signif.stars = getOption("show.signif.stars"), ...)
{
cat("\nCall:\n")
cat(paste(deparse(x$call), sep = "\n", collapse = "\n"),
"\n\n", sep = "")
resid <- x$residuals
df <- x$df
rdf <- df[2]
cat(if (!is.null(x$w) && diff(range(x$w)))
"Weighted ", "Residuals:\n", sep = "")
if (rdf > 5) {
etc.La commande deparse permet d'avoir une chaine de caractères dont l'évaluation reconstruira l'objet initial (le syntaxe est un peu étrange : si on construisait ces objets à la main, on procèderait un peu différemment).
> deparse(r)
[1] "structure(c(-0.0417043206282872, -0.995221281426497, -0.0463357461197109, "
[2] "-0.075155500585547, 0.0448216562696701, 0.058760027857223, -0.97685747990989, "
[3] "-0.200718066450337, 0.0798906594208107, -0.0675697350838044, "
[4] "-0.200546287353865, 0.974080592182492, 0.994921731246978, -0.0389382976351601, "
[5] "0.0581691430589318, -0.0723250196376096), .Dim = c(4, 4), .Dimnames = list("
[6] " c(\"Murder\", \"Assault\", \"UrbanPop\", \"Rape\"), c(\"Comp.1\", \"Comp.2\", "
[7] " \"Comp.3\", \"Comp.4\")), class = \"loadings\")"
> cat(deparse(r)); cat("\n")
structure(c(-0.0417043206282872, -0.995221281426497,
-0.0463357461197109, -0.075155500585547, 0.0448216562696701,
0.058760027857223, -0.97685747990989, -0.200718066450337,
0.0798906594208107, -0.0675697350838044, -0.200546287353865,
0.974080592182492, 0.994921731246978, -0.0389382976351601,
0.0581691430589318, -0.0723250196376096), .Dim = c(4, 4), .Dimnames
= list( c("Murder", "Assault", "UrbanPop", "Rape"), c("Comp.1",
"Comp.2", "Comp.3", "Comp.4")), class = "loadings")
Il y a plusieurs manières de coder les résultats d'une expérience.
Premier exemple : on a mesuré plusieurs variables qualitatives sur plusieurs (centaines de) sujets. On peut représenter ces données sous forme d'un tableau, une ligne par sujet, une colonne par varible. On peut aussi représenter ces données sous forme de tableau de contingence (c'est pertinent s'il y a peu de variables, mais inutilisable sinon : le tableau aura trop de cases nulles).
Comment passer d'une formulation à l'autre ?
Dans un sens, il y a la commande "table", qui calcule un tableau de contingence.
n <- 1000 x1 <- factor( sample(1:3, n, replace=T), levels=1:3 ) x2 <- factor( sample(LETTERS[1:5], n, replace=T), levels=LETTERS[1:5] ) x3 <- factor( sample(c(F,T),n,replace=T), levels=c(F,T) ) d <- data.frame(x1,x2,x3) r <- table(d)
Il vient :
> r , , x3 = FALSE x2 x1 A B C D E 1 27 45 31 38 25 2 41 33 30 35 33 3 33 30 28 35 39 , , x3 = TRUE x2 x1 A B C D E 1 26 30 28 42 29 2 35 33 22 37 40 3 42 31 31 36 35
Il y a aussi une commande ftable, qui présente le résultat différemment.
> ftable(d)
x3 FALSE TRUE
x1 x2
1 A 27 26
B 45 30
C 31 28
D 38 42
E 25 29
2 A 41 35
B 33 33
C 30 22
D 35 37
E 33 40
3 A 33 42
B 30 31
C 28 31
D 35 36
E 39 35Voyons maintenant comment transformer un tableau de contingence en data.frame.
Cas 1 : table uni-dimensionnelle
n <- 100 k <- 10 x <- factor( sample(LETTERS[1:k], n, replace=T), levels=LETTERS[1:k] ) d <- table(x) factor( rep(names(d),d), levels=names(d) )
Cas 2 : table bi-dimensionnelle
n <- 100
k <- 4
x1 <- factor( sample(LETTERS[1:k], n, replace=T), levels=LETTERS[1:k] )
x2 <- factor( sample(c('x','y','z'),n,replace=T), levels=c('x','y','z') )
d <- data.frame(x1,x2)
d <- table(d)
y2 <- rep(colnames(d)[col(d)], d)
y1 <- rep(rownames(d)[row(d)], d)
dd <- data.frame(y1,y2)Cas 3 : cas général
n <- 1000
x1 <- factor( sample(1:3, n, replace=T), levels=1:3 )
x2 <- factor( sample(LETTERS[1:5], n, replace=T), levels=LETTERS[1:5] )
x3 <- factor( sample(c(F,T),n,replace=T), levels=c(F,T) )
d <- data.frame(x1,x2,x3)
r <- table(d)
# Une fonction qui généralise row et col en dimension supérieures.
foo <- function (r,i) {
d <- dim(r)
rep(
rep(1:d[i], each=prod(d[0:(i-1)])),
prod(d[(i+1):(length(d)+1)], na.rm=T)
)
}
k <- length(dimnames(r))
y <- list()
for (i in 1:k) {
y[[i]] <- rep( dimnames(r)[[i]][foo(r,i)], r )
}
d <- data.frame(y)
colnames(d) <- LETTERS[1:k]
# Vérification
r - table(d)
Autre exemple : on fait trois fois la même expérience sur le même individu. On peut choisir de mettre une ligne par sujet, avec donc plusieurs résultats
subject, result1, result2, result3
on peut aussi choisir de mettre une ligne par expérience, avec le numéro du sujet, le numéro de l'expérience et le résultat.
subject, result, retry
Exercice : écrire des fonctions qui effectuent la conversion d'une forme vers l'autre. On pourra utiliser la commande split, qui permet de séparer des données selon la valeur d'un facteur.
Autre exemple : même situation, mais cette fois-ci, le nombre de répétitions par individu n'est pas constant. On peut représenter les données par une liste de vecteurs.
n <- 100 k <- 10 subject <- factor( sample(1:k,n,replace=T), levels=1:k ) x <- rnorm(n) d1 <- data.frame(subject, x) # Data.frame vers liste de vecteurs d2 <- split(d1$x, d1$subject) # Liste de vecteurs vers data.frame rep(names(d2), sapply(d2, length))
A FAIRE : voir s'il ne fait pas fusionner cette section avec la suivante.
En SQL (le langage avec lequel on s'adresse aux bases de données -- pour simplifier, on va considérer qu'une base de donnée, c'est un (ensemble) de data.frame(s)), on veut souvent appliquer une fonction (par exemple la moyenne, la somme) à des groupes d'enregistrements (un enregistrement, c'est une ligne dans le data.frame). Par exemple, si on stocke sa comptabilité personnelle dans une base de données, en spécifiant, pour chaque dépense, le montant et la nature (loyer, nourriture, transport, impôts, livres, cinéma, etc.),
montant nature ------------------ 390 loyer 4.90 cinéma 6.61 nourriture 10.67 nourriture 6.40 livres 14.07 nourriture 73.12 livres 4.90 cinéma
on peut vouloir calculer la somme des dépenses de chaque type. En SQL, on dirait :
A FAIRE
On peut faire la même chose avec R :
A FAIRE ?by
A FAIRE : un exemple avec des données réelles ou simulées.
A FAIRE : mettre ça un peu plus haut ?
La commande apply permet d'appliquer une fonction (par exemple, la moyenne, des quartiles, etc.) à chaque colonne (ou ligne) d'un Data Frame (ou d'une matrice).
> options(digits=4)
> df <- data.frame(x=rnorm(20),y=rnorm(20),z=rnorm(20))
> apply(df,2,mean)
x y z
0.04937 -0.11279 -0.02171
> apply(df,2,range)
x y z
[1,] -1.564 -1.985 -1.721
[2,] 1.496 1.846 1.107Ca marche aussi en dimensions supérieures. Le deuxième argument indique les dimensions qui restent, i.e., les dimensions qu'on utilise pour couper les données en tranches.
> options(digits=2)
> m <- array(rnorm(10^3), dim=c(10,10,10))
> a <- apply(m, 1, mean)
> a
[1] 0.060 -0.027 0.037 0.160 0.054 0.012 -0.039 -0.064 -0.013 0.061
> b <- apply(m, c(1,2), mean)
> b
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] -0.083 -0.7297 0.547 0.283 0.182 -0.409 -0.0029 0.170 -0.131 0.7699
[2,] -0.044 0.3618 -0.206 -0.095 0.062 -0.568 -0.4841 0.334 0.362 0.0056
[3,] 0.255 0.2359 -0.331 0.040 0.213 -0.547 -0.1852 0.492 -0.257 0.4525
[4,] -0.028 0.7422 0.417 -0.088 0.205 -0.521 -0.1981 0.042 0.604 0.4244
[5,] -0.085 0.3461 0.047 0.683 -0.018 -0.173 0.1825 -0.826 -0.037 0.4153
[6,] -0.139 -0.4761 0.276 0.174 0.145 0.232 -0.1194 -0.010 0.176 -0.1414
[7,] -0.139 0.0054 -0.328 -0.264 0.078 0.496 0.2812 -0.336 0.124 -0.3110
[8,] -0.060 0.1291 0.313 -0.199 -0.325 0.338 -0.2703 0.166 -0.133 -0.5998
[9,] 0.091 0.2250 0.155 -0.277 0.075 -0.044 -0.4169 0.050 0.200 -0.1849
[10,] -0.157 -0.3316 -0.103 0.373 -0.034 0.116 0.0660 0.249 -0.040 0.4689
> apply(b, 1, mean)
[1] 0.060 -0.027 0.037 0.160 0.054 0.012 -0.039 -0.064 -0.013 0.061La fonction tapply permet de regroupper les observations selon la valeur d'un (ou plusieurs) facteurs et d'appliquer une fonction (moyenne, etc.) à chacun des groupes ainsi obtenus. La commande by fait à peu près la même chose.
> tapply(1:20, gl(2,10,20), sum) 1 2 55 155 > by(1:20, gl(2,10,20), sum) INDICES: 1 [1] 55 ------------------------------------------------------------ INDICES: 2 [1] 155
La fonction sapply applique une fonction à chaque élément d'une liste (ou d'un vecteur, etc.) et renvoie si possible un vecteur. La fonction lapply fait la même chose, mais renvoie une liste.
> x <- list(a=rnorm(10), b=runif(100), c=rgamma(50,1))
> lapply(x,sd)
$a
[1] 1.041
$b
[1] 0.294
$c
[1] 1.462
> sapply(x,sd)
a b c
1.041 0.294 1.462En particulier, la commande sapply permet d'appliquer une fonction à chaque colonne d'un data.frame sans avoir à préciser le numéro des dimensions à garder comme avec la commande apply.
A FAIRE : la commande split
La commande split permet de découper des données : exactement comme tapply, sauf qu'on n'applique aucune fonction derrière.
A FAIRE : match
?match # Get the 2's and 4's x[as.logical( match(x, c(2,4), nomatch=0) )] A FAIRE ?nchar ?substring A FAIRE
On peut aussi manipuler des dates :
library(help=date) library(help=chron)
Beaucoup de commande R manipulent des vecteurs ou des tableaux, permettant ainsi une programmation (presque) sans boucle, une programmation parallèle. Ainsi, les calculs se font beaucoup plus rapidement que si la boucle était écrite explicitement (car R est un langage interprété). Le style de programmation qui en résulte est très différent de ce à quoi on est habitué : voici donc quelques exercices pour vous rôder. Nous aurons besoin des fonctions de manipulation de tableaux que nous avons vues (en particulier, apply).
Soit x un tableau. Calculer la somme de chacune de ses lignes, la somme de chacune de ses colonnes. Si x est le tableau d'effectifs de deux variables qualitatives, calculer le tableau des effectifs théoriques sous l'hypothèse d'indépendance de ces deux variables. Si vous savez ce que c'est, calculer le Chi^2 correspondant.
# Pour être sûr de ne pas confondre lignes et colonnes, je ne prends # pas un tableau carré. n <- 4 m <- 5 x <- matrix( rpois(n*m,10), nr=n, nc=m ) rownames(x) <- 1:n colnames(x) <- LETTERS[1:m] x apply(x,1,sum) apply(x,2,sum) # Effectifs théoriques y <- matrix(apply(x,1,sum),nr=n,nc=m) * matrix(apply(x,2,sum),nr=n,nc=m,byrow=T) / sum(x) # Effectifs théoriques y <- apply(x,1,sum) %*% t(apply(x,2,sum)) / sum(x) # On peut aussi calculer le Chi2 à la main sum((x-y)^2/y) # Vérification chisq.test(x)$statistic
Soit x un vecteur booléen. Compter le nombre de suites de 0 (par exemple, dans 00101001010110, il y a 6 suites de 0). Compter le nombre de suites de 1. Compter le nombre de suites. Même question avec un facteur à plus de deux valeurs.
n <- 50
x <- sample(0:1, n, replace=T, p=c(.2,.8))
# Nombre de suites
sum(abs(diff(x)))+1
# Nombre de suites de 1
f <- function (x, v=1) { # Si quelqu'un a quelque chose de plus simple...
x <- diff(x==v)
x <- x[x!=0]
if(x[1]==1)
sum(x==1)
else
1+sum(x==1)
}
f(x,1)
# Nombre de suites de 0
f(x,0)
n <- 50
k <- 4
x <- sample(1:k, n, replace=T)
# Avec une boucle
s <- 0
for (i in 1:4) {
s <- s + f(x,i)
}
s
# Sans boucle (peu lisible)
a <- apply(matrix(1:k,nr=1,nc=k), 2, function (i) { f(x,i) } )
a
sum(a)Dans un vecteur binaire de taille n, trouver la position des suites de 1 de longueur supérieure à k.
n <- 100 k <- 10 M <- sample(0:1, n, replace=T, p=c(.2,.8)) x <- c(0,M,0) # Début des suites de 1 deb <- which( diff(x) == 1 ) # Fin des suites de 1 fin <- which( diff(x) == -1 ) -1 # Longueur des suites de 1 long <- fin - deb # Position de début et de fin des suites de 1 de longueur au moins k cbind(deb,fin)[ fin-deb >= k, ]
Exercice : même question, mais on cherche des suites de 1 de longueur au moins k dans les lignes d'une matrice n*m. Présenter le résultat sous forme de tableau.
foo <- function (M,k) {
x <- c(0,M,0)
deb <- which( diff(x) == 1 )
fin <- which( diff(x) == -1 ) -1
cbind(deb,fin)[ fin-deb >= k, ]
}
n <- 50
m <- 50
M <- matrix( sample(0:1, n*m, replace=T, prob=c(.2,.8)), nr=n, nc=m )
res <- apply(M, 1, foo, k=10)
# Ajout du numéro de la ligne (ça n'est pas très élégant : si
# quelqu'un a mieux...
i <- 0
res <- lapply(res, function (x) {
x <- matrix(x, nc=2)
i <<- i+1
#if (length(x)) {
cbind(ligne=rep(i,length(x)/2), deb=x[,1], fin=x[,2])
#} else {
# x
#}
})
# Présentation du résultat sous forme de tableau
do.call('rbind', res) # Mais il manque les numéros de ligne.A FAIRE : vérifier que je parle de la fonction do.call quelque part dans ce document...
A FAIRE : problème de langue...
Let r be the return of a financial asset. The clustered return is the accumulative return for a sequence of returns of the same sign. The trend number is the number of steps in such a sequence. The average return is their ratio. Compute these quantities.
data(EuStockMarkets)
x <- EuStockMarkets
# We aren't interested in the spot prices, but in the returns
# return[i] = ( price[i] - price[i-1] ) / price[i-1]
y <- apply(x, 2, function (x) { diff(x)/x[-length(x)] })
# We normalize the data
z <- apply(y, 2, function (x) { (x-mean(x))/sd(x) })
# A single time series
r <- z[,1]
# The runs
f <- factor(cumsum(abs(diff(sign(r))))/2)
r <- r[-1]
accumulative.return <- tapply(r, f, sum)
trend.number <- table(f)
boxplot(abs(accumulative.return) ~ trend.number, col='pink',
main="Accumulative return")
boxplot(abs(accumulative.return)/trend.number ~ trend.number,
col='pink', main="Average return")
op <- par(mfrow=c(2,2))
for (i in 1:4) {
r <- z[,i]
f <- factor(cumsum(abs(diff(sign(r))))/2)
r <- r[-1]
accumulative.return <- tapply(r, f, sum)
trend.number <- table(f)
boxplot(abs(accumulative.return) ~ trend.number, col='pink')
}
par(op)
op <- par(mfrow=c(2,2))
for (i in 1:4) {
r <- z[,i]
f <- factor(cumsum(abs(diff(sign(r))))/2)
r <- r[-1]
accumulative.return <- tapply(r, f, sum)
trend.number <- table(f)
boxplot(abs(accumulative.return)/trend.number ~ trend.number, col='pink')
}
par(op)
Soit M une matrice n*m (représentant une image en niveaux de gris) ; calculer la valeur moyenne de chaque quadripixel.
data(volcano) M <- volcano n <- dim(M)[1] m <- dim(M)[2] M1 <- M [1:(n-1),] [,1:(m-1)] M2 <- M [2:n,] [,1:(m-1)] M3 <- M [1:(n-1),] [,2:m] M4 <- M [2:n,] [,2:m] # Avec des quadripixels qui se chevauchent M0 <- (M1+M2+M3+M4)/4 # Avec des quadripixels qui ne se chevauchent pas nn <- floor((n-1)/2) mm <- floor((m-1)/2) M00 <- M0 [2*(1:nn),] [,2*(1:mm)] op <- par(mfrow=c(2,2)) image(M, main="Image initiale") image(M0, main="Quadripixels chevauchants") image(M00, main="Quadripixels non chevauchants") par(op)
Calculer une matrice de Van der Monde.
outer(x, 0:n, '^')
Tracer un graphe à l'aide de sa matrice d'incidence
n <- 100
m <- matrix(runif(2*n),nc=2)
library(ape)
r <- mst(dist(m)) # La matrice d'incidence (de l'arbre couvrant
# minimal des points précédents).
plot(m)
n <- dim(r)[1]
w <- which(r!=0)
i <- as.vector(row(r))[w]
j <- as.vector(col(r))[w]
segments( m[i,1], m[i,2], m[j,1], m[j,2], col='red' )
A FAIRE : trouver d'autres exercices...
Les opérateurs suivants ont le sens auquel vous pensez (mais ils ont tendance à manipuler des vecteurs) :
+ * / - ^ < <= > >= == !=
Les opérateurs logiques sont !, & et | (mais on peut écrire && ou || à la place de | et & : le résultat est alors toujours un scalaire, même si les arguments sont des vecteurs).
L'opérateur : (deux-points) permet de créer des vecteurs.
> -5:7 [1] -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7
L'opérateur [ permet de récupérer un (ou plusieurs) éléments d'un tableau.
> x <- floor(10*runif(10)) > x [1] 3 6 5 1 0 6 7 8 5 8 > x[3] [1] 5 > x[1:3] [1] 3 6 5 > x[c(1,2,5)] [1] 3 6 0
L'opérateur $ permet de récupérer un élément d'une liste, sans avoir à mettre son nom entre guillemets, contrairement à l'opérateur [[. L'intéret de l'opérateur [[ est que son argument peut être une variable.
> op <- par() > op$col [1] "black" > op[["col"]] [1] "black" > a <- "col" > op[[a]] [1] "black"
L'assignation se note <-.
x <- 1.17 y <- c(1, 2, 3, 4)
Le produit matriciel se note %*%, le produit tensoriel (aka produit de Kronecker) %x%.
> A <- matrix(c(1,2,3,4), nr=2, nc=2)
> J <- matrix(c(1,0,2,1), nr=2, nc=2)
> A
[,1] [,2]
[1,] 1 3
[2,] 2 4
> J
[,1] [,2]
[1,] 1 2
[2,] 0 1
> J %x% A
[,1] [,2] [,3] [,4]
[1,] 1 3 2 6
[2,] 2 4 4 8
[3,] 0 0 1 3
[4,] 0 0 2 4L'opérateur %o% sert à fabriquer des tables de multiplications (il appelle la fonction outer avec la multiplication).
> A <- 1:5 > B <- 11:15 > names(A) <- A > names(B) <- B > A %o% B 11 12 13 14 15 1 11 12 13 14 15 2 22 24 26 28 30 3 33 36 39 42 45 4 44 48 52 56 60 5 55 60 65 70 75
La division euclidienne est %/%, son reste est %%.
> 1234 %% 3 [1] 1 > 1234 %/% 3 [1] 411 > 411*3 + 1 [1] 1234
L'appartenance à un "ensemble" se note %in%.
> 17 %in% 1:100 [1] TRUE > 17.1 %in% 1:100 [1] FALSE
L'opérateur ~ est utilisé pour décrire des modèles : nous en parlerons plus loin.
Pour plus de détails (et pour les opérateurs que j'ai passé sous silence, comme <<- -> ->> @ :: _ =), voir le manuel :
?"+" ?"<" ?"<-" ?"!" ?"[" ?Syntax ?kronecker ?match library(methods) ?slot
On peut s'amuser à définir ses propres opérateurs : ce sont juste des fonction à deux arguments dont le nom commence et se termine par %. L'exemple suivant est tiré du manuel.
> "%w/o%" <- function(x,y) x[!x %in% y] > (1:10) %w/o% c(3,7,12) [1] 1 2 4 5 6 8 9 10
R est en fait un langage de programmation : à ce titre, on dispose des structures de contrôle habituelles (boucles, conditioinnelles, récursivité, etc.).
Conditionnelle :
if(...) {
...
} else {
...
}On peut utiliser les conditionnelles à l'intérieur d'autres constructions.
x <- if(...) 3.14 else 2.71
On peut aussi construire des vecteurs à l'aide de conditionnelles.
x <- rnorm(100) y <- ifelse(x>0, 1, -1) z <- ifelse(x>0, 1, ifelse(x<0, -1, 0))
Branchement :
x <- letters[floor(1+runif(1,0,4))]
y <- switch(x,
a='Bonjour',
b='Gutten Tag',
c='Hello',
d='Konnichi wa',
)Boucle for :
for (i in 1:10) {
...
if(...) { next }
...
if(...) { break }
...
}Boucle while :
while(...) {
...
}Boucle repeat :
repeat {
...
if(...) { break }
...
}
R se ratache à la famille des langages fonctionnels (Lisp, OCaML, mais aussi Python) : la notion de fonction est centrale (en particulier, si on en a besoin, il est très facile d'écrire des fonctions qui prennent en argument d'autres fonctions).
On définit une fonction ainsi :
f <- function(x) {
x^2 + x + 1
}La valeur de retour est la dernière valeur calculée (mais on peut aussi utiliser la commande return).
On peut donner des valeurs par défaut aux arguments.
f <- function(x, y=3) { ... }Lors de l'appel de la fonction, on peut utiliser le nom des arguments, sans tenir compte de leur ordre.
f(y=4, x=3.14)
A la fin des arguments, on peut mettre trois points de suspension, qui représentent tous les arguments que l'on n'a pas précisés, et que l'on peut passer en argument à d'autres fonctions.
f <- function(x, ...) {
plot(x, ...)
}Les fonctions n'ont pas d'effet de bord : toutes les modifications sont locales. En particulier, on ne peut pas écrire de fonction qui modifie une variable globale. (En fait, si, on peut : voir un peu plus loin la partie "bidouillages".)
Vincent Zoonekynd
<zoonek@math.jussieu.fr>
latest modification on Wed Oct 13 22:32:30 BST 2004