The French version of this document is no longer maintained: be sure to check the more up-to-date English version.
Exemples de données
Données univariées quantitatives
Données univariées ordonnées
Données univariées qualitatives
Données bivariées quantitatives
Données bivariées qualitative/quantitative
Données bivariées qualitatives
Effets graphiques
Données multivariées quantitatives
Données multivariées en partie qualitatives
Treillis
Dans cette partie, nous montrons comment transformer des données en graphiques avec R, qu'il s'agisse de graphiques simples pour des données univariées ou bivariées, ou des graphiques dont la compréhension requiert un peu d'algèbre linéaire ou des algorithmes non triviaux.
Il y a plein d'exemples de données qui viennent avec R et avec lesquelles on peut jouer.
En voici un.
> ?faithful
> data(faithful)
> faithful
eruptions waiting
1 3.600 79
2 1.800 54
3 3.333 74
...
270 4.417 90
271 1.817 46
272 4.467 74
> 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 ...
Pour regarder quelles sont les données présentes dans un paquetage :
data(package='ts')
Data sets in package `ts':
AirPassengers Monthly Airline Passenger Numbers 1949-1960
austres Quarterly Time Series of the Number of
Australian Residents
beavers Body Temperature Series of Two Beavers
BJsales Sales Data with Leading Indicator.
EuStockMarkets Daily Closing Prices of Major European Stock
Indices, 1991-1998.
JohnsonJohnson Quarterly Earnings per Johnson & Johnson Share
LakeHuron Level of Lake Huron 1875--1972
lh Luteinizing Hormone in Blood Samples
lynx Annual Canadian Lynx trappings 1821--1934
Nile Flow of the River Nile
nottem Average Monthly Temperatures at Nottingham,
1920--1939
sunspot Yearly Sunspot Data, 1700--1988. Monthly
Sunspot Data, 1749--1997.
treering Yearly Treering Data, -6000--1979.
UKDriverDeaths Road Casualties in Great Britain 1969--84
UKLungDeaths Monthly Deaths from Lung Diseases in the UK
UKgas UK Quarterly Gas Consumption
USAccDeaths Accidental Deaths in the US 1973--1978
WWWusage Internet Usage per MinuteLes données présentes dans le paquetage "base" :
Data sets in package `base':
Formaldehyde Determination of Formaldehyde concentration
HairEyeColor Hair and eye color of statistics students
InsectSprays Effectiveness of insect sprays
LifeCycleSavings Intercountry life-cycle savings data
OrchardSprays Potency of orchard sprays
PlantGrowth Results from an experiment on plant growth
Titanic Survival of passengers on the Titanic
ToothGrowth The effect of vitamin C on tooth growth in
guinea pigs
UCBAdmissions Student admissions at UC Berkeley
USArrests Violent crime statistics for the USA
USJudgeRatings Lawyers' ratings of state judges in the US
Superior Court
USPersonalExpenditure Personal expenditure data
VADeaths Death rates in Virginia (1940)
airmiles Passenger miles on US airlines 1937-1960
airquality New York Air Quality Measurements
anscombe Anscombe's quartet of regression data
attenu Joiner-Boore Attenuation Data
attitude Chatterjee-Price Attitude Data
cars Speed and Stopping Distances for Cars
chickwts The Effect of Dietary Supplements on Chick
Weights
co2 Moana Loa Atmospheric CO2 Concentrations
discoveries Yearly Numbers of `Important' Discoveries
esoph (O)esophageal Cancer Case-control study
euro Conversion rates of Euro currencies
eurodist Distances between European Cities
faithful Old Faithful Geyser Data
freeny Freeny's Revenue Data
infert Secondary infertility matched case-control
study
iris Edgar Anderson's Iris Data as data.frame
iris3 Edgar Anderson's Iris Data as 3-d array
islands World Landmass Areas
longley Longley's Economic Regression Data
morley Michaelson-Morley Speed of Light Data
mtcars Motor Trend Car Data
nhtemp Yearly Average Temperatures in New Haven CT
phones The Numbers of Telephones
precip Average Precipitation amounts for US Cities
presidents Quarterly Approval Ratings for US Presidents
pressure Vapour Pressure of Mercury as a Function of
Temperature
quakes Earthquake Locations and Magnitudes in the
Tonga Trench
randu Random Numbers produced by RANDU
rivers Lengths of Major Rivers in North America
sleep Student's Sleep Data
stackloss Brownlee's Stack Loss Plant Data
state US State Facts and Figures
sunspots Monthly Mean Relative Sunspot Numbers
swiss Swiss Demographic Data
trees Girth, Height and Volume for Black Cherry
Trees
uspop Populations Recorded by the US Census
volcano Topographic Information on Auckland's Maunga
Whau Volcano
warpbreaks Breaks in Yarn during Weaving
women Heights and Weights of WomenIl y en a aussi beaucoup dans le paquetage MASS (qui correspondent au livre "Modern Applied Statistics with S").
data(package='MASS')
R précise même :
Use `data(package = .packages(all.available = TRUE))' to list the data sets in all *available* packages.
C'est un peu violent, mais on peut prendre tous les exemples et leur faire subir un certain trairement. J'ai parfois utilisé de telles manipulations lors de la rédaction de ces notes pour trouver un exemple vérifiant certaines propriétés.
#!/bin/sh
cd /usr/lib/R/library/
for lib in *(/)
do
if [ -d $lib/data ]
then
(
cd $lib/data
echo 'library('$lib')'
ls
) | grep -vi 00Index | perl -p -e 's#(.*)\..*?$#data($1); doit($1)#'
fi
done > /tmp/ALL.RJe peux avoir une liste de toutes les données comme suit :
doit <- function (...) {}
source("/tmp/ALL.R")
sink("str_data")
ls.str()
q()Je peux chercher des données avec des valeurs extrèmes (outliers) ainsi :
doit <- function (d) {
name <- as.character(substitute(d))
cat(paste("Processing", name, "\n"))
if (!exists(name)) {
cat(paste(" Skipping (does not exist!)", name, "\n"))
} else if (is.vector(d)) {
doit2(d, name)
} else if (is.data.frame(d)) {
for (x in names(d)) {
cat(paste("Processing", name, "/", x, "\n"))
doit2(d[[x]], paste(name,"/",x))
}
} else
cat(paste(" Skipping (unknown reason)", substitute(d), "\n"))
}
doit2 <- function (x,n) {
if( is.numeric(x) )
really.do.it(x,n)
else
cat(" Skipping (non numeric)\n")
}
really.do.it <- function (x,name) {
x <- x[!is.na(x)]
m <- median(x)
i <- IQR(x)
n1 <- sum(x>m+1.5*i)
n2 <- sum(x<m-1.5*i)
n <- length(x)
p1 <- round(100*n1/n, digits=0)
p2 <- round(100*n2/n, digits=0)
if( n1+n2>0 ) {
boxplot(x, main=name)
cat(paste(" OK ", n1+n2, "/", n, " (",p1,"%, ",p2,"%)\n", sep=''))
}
}
source("ALL.R", echo=F)On peut aussi procéder ainsi pour tester du code.
On parle de données univariées quand il n'y a qu'une seule variable, bivariées quand il y en a deux, multivariées quand il y en a plus.
On parle de données quantitatives ou numériques, quand cette variable contient des nombres, avec lesquels faire de l'arithmétique a un sens : par exemple des températures, mais pas des numéros de département.
On parle de données qualitatives en cas contraire : par exemple une variable dont les valeurs sont "oui" ou "non", dont les valeurs sont des couleurs, ou dont les valeurs sont des numéros de département.
On rencontre parfois des données qualitatives ordonnées. Par exemple, une variable dont les valeurs seraient "jamais", "rarement", "de temps à autres", "souvent", "tout le temps". Ces données peuvent provenir de la répartition en classes de données quantitatives.
On peut résumer ces données à l'aide de quelques nombres : la moyenne, le minimum, le maximum, la médianne (on met les nombres dans l'ordre et on prend celui du milieu) les quartiles (idem, mais on prend les nombres à un quard du début ou de la fin).
> summary(faithful$eruptions) Min. 1st Qu. Median Mean 3rd Qu. Max. 1.600 2.163 4.000 3.488 4.454 5.100 > mean(faithful$eruptions) [1] 3.487783
Il faut toujours regarder les données avec un oeil critique : il faut en particulier se demander si les valeurs extrémales sont normales ou si elles ne proviendraient pas d'une quelconque erreur.
On peut aussi regarder diverses mesures de la dispersion, comme l'écart moyen à la médiane (MAD), l'écart-type ("standard deviation", en anglais) ou l'intervalle inter-quartiles (IQR).
> mad(faithful$eruptions) [1] 0.9510879 > sd(faithful$eruptions) [1] 1.141371 > IQR(faithful$eruptions) [1] 2.2915
Je profite de l'occasion pour signaler les liens entre moyenne, variance, médiane et écart moyen à la moyenne. Quand on a un ensemble de nombres x1, x2, ... xn, on peut essayer de trouver le réel m qui minimise la somme
abs(m-x1) + abs(m-x2) + ... + abs(m-xn).
C'est la médiane. On peut aussi essayer de trouver le réel m qui minimise la somme
(m-x1)^2 + (m-x2)^2 + ... + (m-xn)^2.
C'est la moyenne. (Cette propriété de la moyenne s'appelle la propriété des moindres carrés.)
D'où la définition suivante : la variance d'une série statistique x1, x2, ..., xn est la moveyye des carrés des écarts à la moyenne :
(m-x1)^2 + (m-x2)^2 + ... + (m-xn)^2
Var(x) = --------------------------------------.
nAu début, j'avais du mal à comprendre la pertinence de la variance comme mesure de la dispersion ; la moyenne des valeurs absolues des différences avec la moyenne me semblait plus naturelle -- en effet, pourquoi introduire ces carrés ? La définition précédente de la moyenne peut donc servir de motivation de la notion de variance. Il y a d'autres motivations. L'une est que c'est facile à calculer, à la main, de manière incrémentale. Une autre est que ça intervient dans plein de résultats théoriques (inégalité de Bienaimé-Tchebichev, etc.).
Les notions de moyenne et de variance perdent de leur pertinence pour des distributions non symétriques ou pour des distributions avec beaucoup de valeurs extrèmes.
La médiane, l'intervalle interquartile ou l'écart moyen à la médiane sont des mesures robustes du centre ou de la dispersion : ils donnent de bons résultats, même sur un échantillon non normal, contrairement à la moyenne ou à l'écart-type. Par contre, si l'échantillon est normal, ils sont moins efficaces : leur variance est plus grande. Nous reviendrons sur ces différences quand nous parlerons d'estimateurs.
A FAIRE : est-ce le bon endroit (on est juste après la variance, certes, mais on est aussi au milieu de notions simples...)
Nous venons de voir que (pour une série statistique centrée, i.e., de moyenne zéro) la variance est la moyenne des carrés. On peut regarder ce qui se passe avec d'autres exposants : le moment d'ordre k d'une série statistique est la moyenne de ses puissances k. Si on le note M_k, on a
moyenne = M_1 variance = M_2 - M_1^2
Le moment d'ordre 3 s'appelle ********A FAIRE******* (en anglais : skewness). Pour une distribution symétrique (centrée sur zéro), il est nul. Pour vérifier si une distribution est symétrique, et pour quantifier l'écart à la symétrie, il suffit de calculer le troisième moment de la série normalisée.
A FAIRE : exemple.
Le moment d'ordre 4 (kurtosis) permet de voir si la distribution est éloignée d'une distribution normale, et de quantifier cet écart.
On le rencontre par exemple dans l'étude des données financières : on suppose souvent que les données qu'on manipule suivent une loi normale, mais dans ce domaine, ce n'est pas le cas. C'est d'autant plus grave que la différence avec la distribution normale se révèle par des valeurs atypiques anormalement nombreuses. Pour le voir, on calcule la densité de la distribution des retours sur investissement et on la superpose à la distribution d'une loi normale. Pour mieux voir, on a pris une échelle logarithmique. On remarque deux choses : d'une part, la distribution est "plus pointue" que la distribution normale, d'autre part, il y a plus de valeurs atypiques.
op <- par(mfrow=c(2,2))
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) })
for (i in 1:4) {
d <- density(z[,i])
plot(d$x,log(d$y),ylim=c(-5,1),xlim=c(-5,5))
curve(log(dnorm(x)),col='red',add=T)
}
par(op)
On vérifie cela par le calcul :
> apply(z^3,2,mean)
DAX SMI CAC FTSE
-0.4344056 -0.5325112 -0.1059855 0.1651614
> apply(z^4,2,mean)
DAX SMI CAC FTSE
8.579151 8.194598 5.265506 5.751968Alors que pour une distribution normale, on obtiendrait 0 et 3.
> mean(rnorm(100000)^3) [1] -0.003451044 > mean(rnorm(100000)^4) [1] 3.016637
On peut aussi faire plusieurs simulations (comme on vient de le faire) et regarder la distribution des valeurs obtenues : celles qui viennent de nos données sont-elles si grandes que ça ?
Pour le troisième moment, deux des valeurs sont extrèmes, mais les deux autres sont assez normales.
n <- dim(z)[1]
N <- 2000
m <- matrix(rnorm(n*N), nc=N, nr=n)
a <- apply(m^3,2,mean)
b <- apply(z^3,2,mean)
hist(a, col='light blue', xlim=range(c(a,b)),
main="Troisième moment")
points(b, rep(.2*par("usr")[3] + .8*par("usr")[4], length(b)), type='h', col='red',lwd=3)
points(b, rep(.2*par("usr")[3] + .8*par("usr")[4], length(b)), col='red', lwd=3)
Par contre, pour la kurtosis, nos valeurs sont vraiment anormalement élevées.
n <- dim(z)[1]
N <- 2000
m <- matrix(rnorm(n*N), nc=N, nr=n)
a <- apply(m^4,2,mean)
b <- apply(z^4,2,mean)
hist(a, col='light blue', xlim=range(c(a,b)),
main="Distribution attendue de la kurtoris et valeurs observées")
points(b, rep(.2*par("usr")[3] + .8*par("usr")[4], length(b)), type='h', col='red',lwd=3)
points(b, rep(.2*par("usr")[3] + .8*par("usr")[4], length(b)), col='red', lwd=3)
Nous reverrons un peu plus loin comment quantifier cet écart à la normalité -- ce que nous venons de faire s'appelle un "calcul de p-valeur par bootstrap paramétrique".
On peut voir les données qu'on étudie comme une masse informe de nombres, dans laquelle on ne voit rien (c'est pourquoi j'ai précédemment utilisé la commande str, pour n'avoir que les début des données : les inclure en intégralité n'apporterait rien).
Une première manière d'y voir quelque chose consiste à les ordonner : mais on a toujours plus d'une centaine de nombres, on n'y voit rien.
> str( sort(faithful$eruptions) ) num [1:272] 1.60 1.67 1.70 1.73 1.75 ...
Mais on remarque que, dans ces nombres, les deux premiers chiffres sont presque toujours les mêmes. De plus, après ces deux premiers chiffres, il n'en reste qu'un. On peut donc les ranger dans différentes classes, selon les deux premiers chiffres, et indiquer à côté le chiffre restant. C'est ce que les anglo-saxons appellent un "stem-and-leaf plot" (graphique avec une tige et des feuilles). C'est juste une autre manière de représenter le tas de nombres initial.
> stem(faithful$eruptions) The decimal point is 1 digit(s) to the left of the | 16 | 070355555588 18 | 000022233333335577777777888822335777888 20 | 00002223378800035778 22 | 0002335578023578 24 | 00228 26 | 23 28 | 080 30 | 7 32 | 2337 34 | 250077 36 | 0000823577 38 | 2333335582225577 40 | 0000003357788888002233555577778 42 | 03335555778800233333555577778 44 | 02222335557780000000023333357778888 46 | 0000233357700000023578 48 | 00000022335800333 50 | 0370
Et on pourrait même faire ça à la main. Voici un peu plus d'explications sur ce type de représentation :
http://www.shodor.org/interactivate/discussions/steml.html http://davidmlane.com/hyperstat/A28117.html http://www.google.fr/search?q=stem-and-leaf&ie=UTF-8&oe=UTF-8&hl=fr&btnG=Recherche+Google&meta=
Un premier graphique consiste à mettre les points sur un axes, tout simplement.
data(faithful) stripchart(faithful$eruptions)
Néanmoins, s'il y a beaucoup de données, ou s'il y a plusieures observations qui ont la même valeur, le résultat est peu lisible. On peut alors ajouter un peu de bruit, de manière que les points ne se superposent pas.
# Bruit uniquement horizontal stripchart(faithful$eruptions, jitter=TRUE)
stripchart(faithful$eruptions, method='jitter')
C'est un bon exercice de familiarisation des commandes rnorm et autres que de tenter de faire ça soi-même.
my.jittered.stripchart <- function (x) {
x.name <- deparse(substitute(x))
n <- length(x)
plot( runif(n) ~ x, xlab=x.name, ylab='random' )
}
my.jittered.stripchart(faithful$eruptions)
my.jittered.stripchart <- function (x) {
x.name <- deparse(substitute(x))
n <- length(x)
x <- x + diff(range(x))*.05* (-.5+runif(n))
plot( runif(n) ~ x, xlab=paste("jittered", x.name), ylab='random' )
}
my.jittered.stripchart(faithful$eruptions)
On peut aussi représenter les données elles-mêmes, après les avoir triées,
plot(sort(faithful$eruptions))
(les deux paliers horizontaux correspondent aux deux pics de l'histogramme).
En fait, c'est simplement un diagramme de dispersion auquel on a ajouté une dimension.
plot(sort(faithful$eruptions)) rug(faithful$eruptions, side=2)
Cela permet aussi de voir si les données sont discrètes (si on fait un diagramme de dispersion sans ajout de bruit (jitter), on peut ne pas s'en appercevoir).
x <- round( rnorm(100), digits=1 ) plot(sort(x)) rug(jitter(x), side=2)
On peut aussi considérer le graphe des effectifs cumulés (c'est juste le symétrique du précédent).
effectifs.cumules <- function (x) {
x.name <- deparse(substitute(x))
n <- length(x)
plot( 1:n ~ sort(x), xlab=x.name, ylab='effectifs cumulés')
}
effectifs.cumules(faithful$eruptions)
Dans certains cas, les individus portent des noms : on peut les faire figurer sur le graphique (c'est le même type de graphique que précédemment, mais tourné de 90 degrés).
data(islands) dotchart(islands)
dotchart(sort(log(islands)))
D'un point de vue théorique, la courbe des effectifs (ou plutôt des fréquences cumulées) est très importante, mais son interprétation ne saute pas vraiment aux yeux. Dans les exemples suivants, on montre la densité de différentes distributions et les courbes de répartition correspondantes.
op <- par(mfcol=c(2,4))
do.it <- function (x) {
hist(x, probability=T, col='light blue')
lines(density(x), col='red', lwd=3)
x <- sort(x)
q <- ppoints(length(x))
plot(q~x, type='l')
abline(h=c(.25,.5,.75), lty=3, lwd=3, col='blue')
}
n <- 200
do.it(rnorm(n))
do.it(rlnorm(n))
do.it(-rlnorm(n))
do.it(rnorm(n, c(-5,5)))
par(op)
Les boites à moustache sont une version simplifiée de la courbe de répartition : on se contente de mettre les quartiles, i.e., les intersections avec les droites pointillées bleues horizontales ci-dessus.
On peut représenter graphiquement les 5 quartiles (minimum, premier quartile, médiane, troisième quratile et maximum) : on obtient une boite à moustaches.
boxplot(faithful$eruptions, range=0)
On comprend mieux cette appellation si on la représente horizontalement.
boxplot(faithful$eruptions, range=0, horizontal=TRUE)
Sur cet exemple, on voit clairement que les données ne sont pas symétriques : on sait donc que ce serait une mauvaise idée de leur appliquer des procédures statistiques qui supposent que les données sont symétriques, ou même carrément normales.
C'est l'une des utilités de ce type de graphique : regarder si les données sont symétriques.
Cette représentation des quartiles est plus simple et plus directement compréhensible que la suivante, en termes d'aires :
op <- par(mfrow=c(1,2))
do.it <- function (x) {
d <- density(x)
plot(d, type='l')
q <- quantile(x)
do.it <- function (i, col) {
x <- d$x[i]
y <- d$y[i]
polygon( c(x,rev(x)), c(rep(0,length(x)),rev(y)), border=NA, col=col )
}
do.it(d$x <= q[2], 'red')
do.it(q[2] <= d$x & d$x <= q[3], 'green')
do.it(q[3] <= d$x & d$x <= q[4], 'blue')
do.it(d$x >= q[4], 'yellow')
lines(d, lwd=3)
}
do.it( rnorm(2000) )
do.it( rlnorm(200) )
par(op)
(Dans cet exemple, les quatre aires sont les mêmes : on affirme souvent que l'oeil humain n'est pas capable de comparer des aires, en voici une preuve éclatante.)
Si on ne précise pas l'option << range=0 >>, le graphe va mettre en évidence les points atypiques, i.e., très éloignés ("outliers", en anglais), qui sont au delà de 1.5 fois l'intervalle interquartile. Sur cet exemple, il n'y en a pas.
boxplot(faithful$eruptions, horizontal=TRUE)
Dans certains cas, c'est tout à fait normal.
# Il y en a, c'est génant, mais pas pathologique boxplot(rnorm(500), horizontal=TRUE)
S'il n'y en a que quelques uns, très isolés, il peut s'agit d'erreurs.
x <- c(rnorm(30),20) x <- sample(x, length(x)) boxplot( x, horizontal=T )
library(boot) data(aml) boxplot( aml$time, horizontal=T )
Ils peuvent aussi signaler que la distribution n'est pas du tout normale.
data(attenu) boxplot(attenu$dist, horizontal=T)
On tente alors de se ramener à des données plus normales en appliquant une fonction bien choisie (nous reviendrons sur cela plus loin).
data(attenu) boxplot(log(attenu$dist), horizontal=T)
La présence de données extrèmes est assez génante, car beaucoup de traitements statistiques (le calcul de la moyenne, de l'écart-type, la régression) y sont très sensibles.
> x <- aml$time > summary(x) Min. 1st Qu. Median Mean 3rd Qu. Max. 5.00 12.50 23.00 29.48 33.50 161.00 > y <- x[x<160] > summary(y) Min. 1st Qu. Median Mean 3rd Qu. Max. 5.00 12.25 23.00 23.50 32.50 48.00
C'est l'autre utilité des boites à moustaches : repérer les valeurs extrèmes, qui peuvent être normales (mais dont il faut être conscient, car elles risquent de fausser les calculs ultérieurs) ; qui peuvent résulter d'erreurs, qu'il convient de corriger ; ou qui signale que la distribution n'est pas normale et contient beaucoup de valeurs extrèmes.
En fait, plus l'échantillon est gros, plus on a de valeurs extrèmes.
y <- c(rnorm(10+100+1000+10000+100000)) x <- c(rep(1,10), rep(2,100), rep(3,1000), rep(4,10000), rep(5,100000)) x <- factor(x) plot(y~x, horizontal=T, las=1)
On pourrait dessiner des boites dont les moustaches s'étendent plus loin pour des échantillons plus gros, mais attention : si la présence de valeurs extrèmes dans les gros échantillons est normale, il peuvent toujours avoir un effet de levier non négligeable.
Exercice : dessiner des boites à moustaches dont la longueur des moustaches s'étend plus loin que 1.5 fois l'intervalle interquartiles pour des échantillons plus gros.
On peut décorer cette boite à moustaches d'un intervalle de confiance pour la médiane.
boxplot(faithful$eruptions, notch=TRUE)
data(breslow) boxplot(breslow$n, notch=TRUE)
On peut aussi ajouter un diagramme de dispersion.
boxplot(faithful$eruptions, horizontal=T) rug(faithful$eruption, ticksize=.2)
On peut aussi représenter ces données par un histogramme : on répartit les données dans différentes classes (l'ordinateur peut s'en charger) et, pour chaque classe, on trace une barre verticale d'autant plus haute que la classe comporte d'éléments.
hist(faithful$eruptions)
Ce qui est assez génant, c'est que des choix de classes différents peuvent donner des histogrammes complètement différents.
Il y a d'une part l'effet de la largeur des classes.
hist(faithful$eruptions, breaks=20, col="light blue")
Mais la position des classes peut complètement changer l'histogramme et le faire apparaitre tantôt symétrique, tantôt non, quand bien-même la largeur des classes ne changerait pas. Ainsi, dans les deux histogrammes suivants, le deuxième pic n'a pas l'air symétrique, mais dissymétrie n'est pas dans le même sens...
op <- par(mfrow=c(2,1)) hist(faithful$eruptions, breaks=seq(1,6,.5), col='light blue') hist(faithful$eruptions, breaks=.25+seq(1,6,.5), col='light blue') par(op)
On peut remplacer l'histogramme par le graphe d'une fonction. Si on voit les données comme une somme de masses de Dirac, il suffit de faire la convolution avec un "noyau" bien choisi, par exemple, la densité de la loi gaussienne.
hist(faithful$eruptions,
probability=TRUE, breaks=20, col="light blue")
points(density(faithful$eruptions, bw=.1), type='l', col='red', lwd=3)
Les courbes de densité présentent aussi le premier problème des histogrammes : un choix de noyau différent donne une courbe différente -- mais le second problème disparait.
hist(faithful$eruptions,
probability=TRUE, breaks=20, col="light blue")
points(density(faithful$eruptions, bw=1), type='l', lwd=3, col='black')
points(density(faithful$eruptions, bw=.5), type='l', lwd=3, col='blue')
points(density(faithful$eruptions, bw=.3), type='l', lwd=3, col='green')
points(density(faithful$eruptions, bw=.1), type='l', lwd=3, col='red')
On peut décorer un histogramme de bien des manières. On peut par exemple ajouter un diagramme de dispersion en dessous ou superposer à sa densité une densité normale.
hist(faithful$eruptions,
probability=TRUE, breaks=20, col="light blue")
rug(faithful$eruptions)
points(density(faithful$eruptions, bw=.1), type='l', lwd=3, col='red')
f <- function(x) {
dnorm(x,
mean=mean(faithful$eruptions),
sd=sd(faithful$eruptions),
)
}
curve(f, add=T, col="red", lwd=3, lty=2)
Quand on cherche à voir si des données sont symétriques, on peut mettre les données dans l'ordre et tenter d'apparier la première avec la dernière, plus la seconde avec l'avant-dernière, etc. C'est ce que fait le graphique suivant, qui représente la distance à la médiane du (n-i)-ième point en fonction de celle du i-ième point.
symetry.plot <- function (x,...) {
x <- x[ !is.na(x) ]
x <- sort(x)
x <- abs(x - median(x))
n <- length(x)
nn <- ceiling(n/2)
plot( x[n:(n-nn+1)] ~ x[1:nn] ,
xlab='Distance below median',
ylab='Distance above median',
...)
}
symetry.plot(faithful$eruptions)
abline(0,1)
Le problème, c'est qu'il est assez délicat de voir quand on est << loin >> de la droite : plus l'échantillon contient de points, plus le graphe est proche de la droite. On voit ici qu'avec une centaine de points, on peut quand même en être assez loin.
n <- 100
N <- 1e3
symetry.plot( rnorm(n), pch='.', col='red',
main=paste("Echantillon de", n, "individus") )
for (i in 1:N) {
x <- sort(rnorm(n))
x <- abs(x - median(x))
l <- length(x)
ll <- ceiling(n/2)
points( x[l:(l-ll+1)] ~ x[1:ll], pch='.', col='red' )
}
abline(0,1)
Avec dix points, c'est encore pire...
(C'est probablement pour cette raison que ce genre de graphique est peu utilisé...)
n <- 10
N <- 1e3
symetry.plot( rnorm(n), pch=15, col='red', cex=.5,
main=paste("Echantillon de", n, "individus") )
for (i in 1:N) {
x <- sort(rnorm(n))
x <- abs(x - median(x))
l <- length(x)
ll <- ceiling(n/2)
points( x[l:(l-ll+1)] ~ x[1:ll], pch=15, col='red', cex=.5 )
}
abline(0,1)
Une manière simple de voir si on est << loin >> ou pas consiste donc à superposer le diagramme de symétrie au graphique précédent (en se souvenant qu'il a été construit avec des échantillons normaux, et qu'avec des échantillons symétriques non normaux, les choses peuvent être différentes).
symetry.plot <- function (x, col2='red',...) {
x <- x[ !is.na(x) ]
x <- sort(x)
x <- abs(x - median(x))
n <- length(x)
nn <- ceiling(n/2)
plot( x[n:(n-nn+1)] ~ x[1:nn] ,
xlab='Distance below median',
ylab='Distance above median',
...)
m <- mean(x)
s <- sd(x)
N <- 1000
for (i in 1:N) {
y <- sort(rnorm(n, mean=m, sd=s))
y <- abs(y - median(y))
n <- length(y)
nn <- ceiling(n/2)
points( y[n:(n-nn+1)] ~ y[1:nn], pch='.', col=col2 )
}
points( x[n:(n-nn+1)] ~ x[1:nn], ...)
abline(0,1, ...)
}
symetry.plot(faithful$eruptions)
symetry.plot(rnorm(100))
symetry.plot(runif(100))
Nous venons de présenter une manière (graphique) de voir si une distribution est symétrique. Une autre chose que l'on aime bien, dans les données qu'on est amené à étudier, ce sont les variables normales.
Dans certains cas, il est flagrant que la variable n'est pas normale : c'est le cas de la durée des eruptions du Old Faithful Geyser (il y a deux pics). Dans d'autres cas, c'est moins flagrant. Une première manière de s'en appercevoir consiste à comparer la densité de nos données avec une densité normale.
data(airquality)
x <- airquality[,4]
hist(x, probability=TRUE, breaks=20, col="light blue")
rug(jitter(x, 5))
points(density(x), type='l', lwd=3, col='red')
f <- function(t) {
dnorm(t, mean=mean(x), sd=sd(x) )
}
curve(f, add=T, col="red", lwd=3, lty=2)
On peut aussi voir graphiquement si une variable est normale : il suffit de représenter les quantiles d'une distribution normale en fonction des quantiles de l'échantillon.
Il y a déjà une fonction qui fait cela. (La fonction qqline trace une droite passant par les premier et troisième quartiles.) Dans cet exemple, les données sont à peu près normales, mais on voit qu'elles sont discrètes.
x <- airquality[,4] qqnorm(x) qqline(x)
Voici ce que l'on peut obtenir avec un échantillon normal.
y <- rnorm(100) qqnorm(y) qqline(y)
Maintenant, avec un échantillon non normal.
y <- rnorm(100)^2 qqnorm(y) qqline(y)
Comme précédemment, on peut juxtaposer plusieurs qqplots normaux à notre graphique, pour voir à quel point les données sont éloignées de la normale.
my.qqnorm <- function (x, N=1000, ...) {
op <- par()
x <- x[!is.na(x)]
n <- length(x)
m <- mean(x)
s <- sd(x)
print("a")
qqnorm(x, axes=F, ...)
for (i in 1:N) {
par(new=T)
qqnorm(rnorm(n, mean=m, sd=s), col='red', pch='.',
axes=F, xlab='', ylab='', main='')
}
par(new=T)
qqnorm(x, ...)
qqline(x, col='blue', lwd=3)
par(op)
}
my.qqnorm(rnorm(100))
my.qqnorm(runif(100), main='loi uniforme')
my.qqnorm(exp(rnorm(100)), main='loi log-normale')
my.qqnorm(c(rnorm(50), 5+rnorm(50)), main='deux pics')
my.qqnorm(c(rnorm(50), 20+rnorm(50)), main='deux pics éloignés')
x <- rnorm(100) x <- x + x^3 my.qqnorm(x, main='distribution très étalée')
Voici d'autres exemples de QQplot.
Deux distributions décalées à gauche.
y <- exp(rnorm(100)) qqnorm(y, main='(1) Loi log-normale'); qqline(y, col='red')
y <- rnorm(100)^2 qqnorm(y, ylim=c(-2,2), main="(2) Carré d'une va suivant une loi normale") qqline(y, col='red')
Une distribution décalée à droite.
y <- -exp(rnorm(100)) qqnorm(y, ylim=c(-2,2), main="(3) Opposé d'une v.a. suivant une loi log-normale") qqline(y, col='red')
Une distribution moins dispersée que la distribution normale (ça s'appelle une distribution leptokurtique).
y <- runif(100, min=-1, max=1) qqnorm(y, ylim=c(-2,2), main='(4) Loi uniforme'); qqline(y, col='red')
Une distribution plus dispersée que la distribution normale (ça s'appelle une distribution platykurtique).
y <- rnorm(10000)^3 qqnorm(y, ylim=c(-2,2), main="(5) Cube d'une v.a. suivant une loi normale") qqline(y, col='red')
Une distribution avec plusieurs pics.
y <- c(rnorm(50), 5+rnorm(50)) qqnorm(y, main='(6) Deux pics') qqline(y, col='red')
y <- c(rnorm(50), 20+rnorm(50)) qqnorm(y, main='(7) Deux pics plus éloignés') qqline(y, col='red')
y <- sample(seq(0,1,.1), 100, replace=T) qqnorm(y, main='(7) Loi discrète') qqline(y, col='red')
On peut lire ces graphiques de la manière suivante.
a. Si la loi est plus concentrée qu'une loi normale à gauche de la moyenne, la partie gauche du graphe se trouve au dessus de la droite (exemples 1, 2 et 4 ci-dessus).
b. Si la loi est moins concentrée qu'une loi normale à gauche de la moyenne, la partie gauche du graphe se trouve en dessous de la droite (exemple 3 ci-dessus).
c. Si la loi est plus concentrée qu'une loi normale à droite de la moyenne, la partie droite du graphe se trouve en dessous de la droite (exemples 3 et 4 ci-dessus).
d. Si la loi est moins concentrée qu'une loi normale à droite de la moyenne, la partie droite du graphe se trouve au dessus de la droite (exemples 1 et 2 ci-dessus).
L'exemple 5 ci-dessus s'interprète donc ainsi : la distribution est symétrique ; à gauche de 0, près de 0, elle est plus concentrée qu'une loi normale ; à gauche de 0, loin de 0, elle est moins concentrée qu'une loi normale ; à droite de 0, c'est pareil.
x <- seq(from=0, to=2, length=100)
y <- exp(x)-1
plot( y~x, type='l', xlim=c(-2,2), ylim=c(-2,2), col='red')
lines( x~y, type='l', col='green')
x <- -x
y <- -y
lines( y~x, type='l', col='blue', )
lines( x~y, type='l', col='cyan')
abline(0,1)
legend( -2, 2,
c( "loi moins concentrée que la normale à droite",
"loi plus concentrée que la normale à droite",
"loi moins concentrée que la normale à gauche",
"loi plus concentrée que la normale à gauche"
),
lwd=3,
col=c("red", "green", "blue", "cyan")
)
title(main="Lecture d'un qqplot")
e. Si la distribution est << plus décentrée vers la gauche qu'avec une loi normale >> (plus précisément, si la médiane est inférieure à la moyenne des premier et troisième quartiles) , alors la courbe passe en dessous de la droite au centre du graphique (exemples 1 et 2 ci-dessus).
f. Si << plus décentrée vers la droite qu'avec une loi normale >> (plus précisément, si la médiane est supérieure à la moyenne des premier et troisième quartiles) , alors la courbe passe au dessus de la droite au centre du graphique (exemple 3 ci-dessus).
g. Si la distribution est symétrique (plus précisément : si la médiane est au milieu des premier et troisème quartiles), alors la courbe coupe la droite au centre du dessin (exemples 4 et 5 ci-dessus).
op <- par()
layout( matrix( c(2,2,1,1), 2, 2, byrow=T ),
c(1,1), c(1,6),
)
# Le graphique
n <- 100
y <- rnorm(n)
x <- qnorm(ppoints(n))[order(order(y))]
par(mar=c(5.1,4.1,0,2.1))
plot( y~x, col="blue" )
y1 <- scale( rnorm(n)^2 )
x <- qnorm(ppoints(n))[order(order(y1))]
lines(y1~x, type="p", col="red")
y2 <- scale( -rnorm(n)^2 )
x <- qnorm(ppoints(n))[order(order(y2))]
lines(y2~x, type="p", col="green")
abline(0,1)
# La légende
# Quelqu'un sait-il comment superposer des graphiques ?
# Je parviens uniquement à les juxtaposer...
# (normalement, c'est le paramètre graphique "fig", mais ça ne marche pas...)
# Si : voir un peu plus loin (il faut préciser par(new=T)).
par(bty='n', ann=F)
g <- seq(0,1, length=10)
e <- g^2
f <- sqrt(g)
h <- c( rep(1,length(e)), rep(2,length(f)), rep(3,length(g)) )
par(mar=c(0,4.1,0,0))
boxplot( c(e,f,g) ~ h, horizontal=T,
border=c("red", "green", "blue"),
xaxt='n',
yaxt='n',
)
title(main="Lecture d'un qqplot")
par(op)
On peut réaliser un QQplot à la main, en reprenant sa définition donnée plus haut.
y <- rnorm(100)^2 y <- scale(x) y <- sort(x) x <- qnorm( seq(0,1,length=length(y)) ) plot(y~x) abline(0,1)
Regardons comment la fonction qqnorm est programmée.
> qqnorm
function (y, ...)
UseMethod("qqnorm")
> help.search("qqnorm")
Help files with alias or title matching `qqnorm',
type `help(FOO, package = PKG)' to inspect entry `FOO(PKG) TITLE':
qqnorm(base) Quantile-Quantile Plots
qqnorm.gls(nlme) Normal Plot of Residuals from a gls Object
qqnorm.lme(nlme) Normal Plot of Residuals or Random Effects
from an lme Object
> apropos("qqnorm")
[1] "my.qqnorm" "qqnorm" "qqnorm.default"
> qqnorm.default
function (y, ylim, main = "Normal Q-Q Plot", xlab = "Theoretical Quantiles",
ylab = "Sample Quantiles", plot.it = TRUE, ...)
{
y <- y[!is.na(y)]
if (0 == (n <- length(y)))
stop("y is empty")
if (missing(ylim))
ylim <- range(y)
x <- qnorm(ppoints(n))[order(order(y))]
if (plot.it)
plot(x, y, main = main, xlab = xlab, ylab = ylab, ylim = ylim,
...)
invisible(list(x = x, y = y))
}On peut faire la même chose pour comparer avec d'autres distributions.
qq <- function (y, ylim, quantiles=qnorm,
main = "Q-Q Plot", xlab = "Theoretical Quantiles",
ylab = "Sample Quantiles", plot.it = TRUE, ...)
{
y <- y[!is.na(y)]
if (0 == (n <- length(y)))
stop("y is empty")
if (missing(ylim))
ylim <- range(y)
x <- quantiles(ppoints(n))[order(order(y))]
if (plot.it)
plot(x, y, main = main, xlab = xlab, ylab = ylab, ylim = ylim,
...)
# From qqline
y <- quantile(y, c(0.25, 0.75))
x <- quantiles(c(0.25, 0.75))
slope <- diff(y)/diff(x)
int <- y[1] - slope * x[1]
abline(int, slope, ...)
invisible(list(x = x, y = y))
}
y <- runif(100)
qq(y, quantiles=qunif)
(Les différentes interprétations du qqplot restent valables, mais les points e, f et g ne renseignent plus sur la symétrie de la distribution, mais comparent cette symétrie avec celle de la loi de référence, qui n'est pas nécessairement symétrique.)
Même si ça n'est pas du tout justifié, on utilise parfois un graphique quantile-quantile pour comparer une variable positive avec une demi-gaussienne (nous le verrons un peu plus loin pour examiner les distances de Cook).
A FAIRE : trouver le nom français...
Pour mieux voir les écarts à la normalité dans un graphique quantiles-quantiles, on peut le tourner, de sorte que la droite soit horizontale.
Attention : ce n'est pas simplement une rotation de 45 degrés car on compare généralement avec une distribution normale standard. Cette droite passe par les premier et troisième quartiles.
two.point.line <- function (x1,y1,x2,y2, ...) {
a1 <- (y2-y1)/(x2-x1)
a0 <- y1 - a1 * x1
abline(a0,a1,...)
}
trended.probability.plot <- function (x, q=qnorm) {
n <- length(x)
plot( sort(x) ~ q(ppoints(n)), , xlab='theoretical quantiles', ylab='sample quantiles')
two.point.line(q(.25), quantile(x,.25), q(.75), quantile(x,.75), col='red')
}
detrended.probability.plot <- function (x, q=qnorm) {
n <- length(x)
x <- sort(x)
x1 <- q(.25)
y1 <- quantile(x,.25)
x2 <- q(.75)
y2 <- quantile(x,.75)
a1 <- (y2-y1)/(x2-x1)
a0 <- y1 - a1 * x1
u <- q(ppoints(n))
x <- x - (a0 + a1 * u)
plot(x ~ u, xlab='theoretical quantiles', ylab='sample quantiles')
abline(h=0, col='red')
}
op <- par(mfrow=c(3,2))
x <- runif(20)
trended.probability.plot(x)
detrended.probability.plot(x)
x <- runif(500)
trended.probability.plot(x)
detrended.probability.plot(x)
trended.probability.plot(x, qunif)
detrended.probability.plot(x,qunif)
par(op)
Cette année, j'ai dû suivre un cours d'analyse de données (tests statistiques, régression, analyse en composantes principales, analyse des correspondances, analyse discriminante, analyse hiérarchique -- mais comme j'avais déjà commencé à rédiger ce document, je n'y ai pas appris grand-chose) au cours duquel on nous a présenté (de manière assez confuse) la courbe de concentration de Gini (ou de Lorenz).
C'est cette courbe qui résume des informations de la forme << 20% des malades bénéficient de 90% des dépenses de la sécurité sociale >>. Plaçons-nous par exemple dans la situation suivante
Pourcentage des malades | Pourcentage des dépenses
----------------------------------------------------
20% | 90%
30% | 95%
50% | 99%il s'agit de pourcentages cummulés :
les 20% les plus dépensiers sont responsables de 90% des dépenses les 30% les plus dépensiers sont responsables de 95% des dépenses les 50% les plus dépensiers sont responsables de 99% des dépenses
On peut représenter graphiquement ces informations.
xy <- matrix(c( 0, 0,
.2, .9,
.3, .95,
.5, .99,
1, 1), byrow=T, nc=2)
plot(xy, type='b', pch=15)
polygon(xy, border=F, col='pink')
lines(xy, type='b', pch=15)
abline(0,1,lty=2)
L'aire entre la courbe et la diagonale est d'autant plus grande que la situation est inégalitaire : c'est l'indice de Gini (en fait, pour que l'indice de Gini varie entre 0 (situation égalitaire) et 1 (situation vraiment très inégalitaire), on multiple cette aire par deux).
Voici un autre exemple :
Dans une cellule,
20 gènes sont exprimés 10^5 fois ("exprimé", ça veut dire transcrit en ARNm)
2700 gènes sont exprimés 10^2 fois
280 gènes sont exprimés 10 fois.Convertissons cela en effectifs cumulés :
20 gènes 100000 ARNm 2720 gènes 100100 ARNm 3000 gènes 100110 ARNm
puis en fréquences :
Gènes | ARNm ------------------- 0.7% | 0.9989012 90.7% | 0.9999001 100% | 1
Voici la courbe (la situation est beaucoup plus inégalitaire que la précédente !) :
x <- c(0,20,2720,3000)/3000 y <- c(0,100000,100100,100110)/100110 plot(x,y, type='b', pch=15) polygon(x,y, border=F, col='pink') lines(x,y, type='b', pch=15) abline(0,1,lty=2)
L'exemple classique pour la courbe de Gini est l'étude de l'inégalité des revenus.
Des fonction pour le dessin de la courbe et le calcul de l'indice de Gini se trouvente dans la bibliothèque ineq. (Les courbes obtenues sont en fait les symétriques des précédentes par rapport à la diagonale : ça ne change rien).
La fonction Gini, dans la bibliothèque ineq, calcule l'indice de concentration de Gini. Il n'est bien défini que si la variable statistique étudiée est positive (dans les exemples précédents, cette variable était << les dépenses de sécurité sociale d'un assuré >>, << le nombre de transcrits d'un gène >>, << le salaire d'un individu >> -- nous ne donnions pas directement cette variable, mais simplement les effectifs (ou effectifs cumulés) de sa répartition en classes).
> n <- 500 > library(ineq) > Gini(runif(n)) [1] 0.3241409 > Gini(runif(n,0,10)) [1] 0.3459194 > Gini(runif(n,10,11)) [1] 0.01629126 > Gini(rlnorm(n)) [1] 0.5035944 > Gini(rlnorm(n,0,2)) [1] 0.8577991 > Gini(exp(rcauchy(n,1))) [1] 0.998 > Gini(rpois(n,1)) [1] 0.5130061 > Gini(rpois(n,10)) [1] 0.1702435 library(ineq) op <- par(mfrow=c(2,3)) plot(Lc(runif(n,0,1)), main="uniforme [0,1]", col='red') plot(Lc(runif(n,0,10)), main="uniforme [0,10]", col='red') plot(Lc(runif(n,10,11)), main="uniforme [10,11]", col='red') plot(Lc(rlnorm(n)), main="log-normale", col='red') plot(Lc(rlnorm(n,0,4)), main="log-normale, plus large", col='red') plot(Lc(rcauchy(n,1)), main="Cauchy", col='red') plot(Lc(exp(rcauchy(n,1))), main="log-Cauchy", col='red') plot(Lc(rpois(n,1)), main="Poisson de moyenne 1", col='red') plot(Lc(rpois(n,10)), main="Poisson de moyenne 10", col='red') par(op)
On peut faire le calcul de l'indice de Gini nous-même.
n <- 500 x <- qlnorm((1:(n-1))/n, 1, 2.25) x <- sort(x) i <- (1:n)/n plot(cumsum(x)/sum(x) ~ i, lwd=3, col='red') abline(0,1) 2*sum(i-cumsum(x)/sum(x))/n
Exercice : dans les exemples précédents, les données étaient présentées de deux manières différentes : des effectifs cumulés et des sommes (premiers exemples) ou des réalisations d'une variables aléatoire. Ecrire le code nécessaire pour passer d'une représentation à l'autre.
On ne s'intéresse pas à la valeur numérique, mais au classement. On peut traiter ces données tantôt comme des données quantitatives (d'ailleurs, quand on a une variable quantitative, on peut ranger ces valeurs dans des classes et les remplacer par leur numéro d'ordre), après une éventuelle transformation (de sorte que les données soit distribuées normalement et pas uniformément), tantôt comme des données qualitatives (car, contrairement aux variables quantitatives, elles ne peuvent prendre qu'un nombre fini de valeurs).
Les facteurs ordonnés se construisent exactement comme les facteurs, la commande factor étant simplement remplacée par la commande ordered.
> l <- c('petit', 'moyen', 'grand')
> ordered( sample(l,10,replace=T), levels=l)
[1] grand moyen moyen grand moyen petit moyen petit moyen moyen
Levels: petit < moyen < grandLes graphiques sont les mêmes qu'avec des facteurs, à la seule différence que l'ordre des facteurs n'est pas arbitraire.
# (Graphique trompeur : l'origine n'est pas sur le dessin) data(esoph) dotchart(table(esoph$agegp))
barplot(table(esoph$agegp))
hist(as.numeric(esoph$agegp),
breaks=seq(.5,.5+length(levels(esoph$agegp)),step=1),
col='light blue')
boxplot(as.numeric(esoph$agegp), horizontal=T)
stripchart(jitter(as.numeric(esoph$agegp),2), method='jitter')
plot(table(esoph$agegp), type='b', pch=7)
On a une liste de données non numériques (en fait, ces données peuvent être codées par des nombres, mais une différence de ces nombres n'a pas de sens : par exemple, des numéros de département ou des numéros de réponse à des questions).
Il y a deux manières de présenter ces données : soit un vecteur de chaines de caractères (ou plus précisément de facteurs : quand on l'imprime, on dirait que ce sont des chaines de caractères, mais comme on s'attend à n'avoir qu'un petit nombre de valeurs, souvent répétées, c'est codé de manière un peu plus efficace), soit une table d'effectifs.
> p <- factor(c("oui", "non"))
> x <- sample(p, 100, replace=T)
> str(x)
Factor w/ 2 levels "non","oui": 1 2 2 2 2 1 2 2 2 2 ...
> table(x)
x
non oui
48 52La deuxième représentation est beaucoup plus compacte. Si on n'a que des variables quantitatives, c'est probablement celle qu'on choisira.
> 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 8Par contre, si on a un mélange de variables qualitatives et quantitatives, on préférera la première représentation.
> data(chickwts) > str(chickwts) `data.frame': 71 obs. of 2 variables: $ weight: num 179 160 136 227 217 168 108 124 143 140 ... $ feed : Factor w/ 6 levels "casein","horseb..",..: 2 2 2 2 2 2 2 2 2 2 ...
Mais n'anticipons pas trop : nous allons commencer par nous limiter au cas de données univariées.
Nous avons déjà vu que la commande table permettait d'obtenir la représentation compacte.
> table(chickwts$feed)
casein horsebean linseed meatmeal soybean sunflower
12 10 12 11 14 12Dans l'autre sens, on peut par exemple utiliser la commande rep.
> x <- table(chickwts$feed) > rep(names(x),x) [1] "casein" "casein" "casein" "casein" "casein" "casein" [7] "casein" "casein" "casein" "casein" "casein" "casein" [13] "horsebean" "horsebean" "horsebean" "horsebean" "horsebean" "horsebean" [19] "horsebean" "horsebean" "horsebean" "horsebean" "linseed" "linseed" [25] "linseed" "linseed" "linseed" "linseed" "linseed" "linseed" [31] "linseed" "linseed" "linseed" "linseed" "meatmeal" "meatmeal" [37] "meatmeal" "meatmeal" "meatmeal" "meatmeal" "meatmeal" "meatmeal" [43] "meatmeal" "meatmeal" "meatmeal" "soybean" "soybean" "soybean" [49] "soybean" "soybean" "soybean" "soybean" "soybean" "soybean" [55] "soybean" "soybean" "soybean" "soybean" "soybean" "sunflower" [61] "sunflower" "sunflower" "sunflower" "sunflower" "sunflower" "sunflower" [67] "sunflower" "sunflower" "sunflower" "sunflower" "sunflower"
data(HairEyeColor) x <- apply(HairEyeColor, 2, sum) barplot(x) title(main="Graphique en colonnes")
barplot(x, col=1, density=c(3,7,11,20), angle=c(45,-45,45,-45)) title(main="Graphique en colonnes")
x <- apply(HairEyeColor, 2, sum)
barplot(as.matrix(x), legend.text=TRUE)
title("Graphique en barre")
barplot(as.matrix(x), horiz=TRUE, col=rainbow(length(x)), legend.text=TRUE) title(main="Graphique en barre")
C'est simplement un graphique en colonnes, dans lequel les données sont mises dans l'ordre (l'ordre des effectifs décroissants). S'il n'y a que très peu de valeurs possibles, ça n'apporte pas grand-chose (d'autant plus que dans les exemples précédents, les données étaient déjà dans l'ordre : nous faisions des disgrammes de Paréto sans le savoir...).
data(attenu) op <- par(las=2) # Écrire les labels perpendiculairement aux axes barplot(table(attenu$event)) title(main="Graphique en colonnes") par(op)
op <- par(las=2) barplot(rev(sort(table(attenu$event)))) title(main="Graphique de Paréto") par(op)
Souvent, on rajoute les effectifs cumulés.
# Avec la commande barplot, je n'y arrive pas.
pareto <- function (x) {
op <- par(mar = c(5, 4, 4, 5) + 0.1)
par(las=2)
if( ! inherits(x, "table") ) {
x <- table(x)
}
x <- rev(sort(x))
plot( x, type='h', axes=F, lwd=16 )
axis(2)
points( x, type='h', lwd=12, col=heat.colors(length(x)) )
y <- cumsum(x)/sum(x)
par(new=T)
plot(y, type="b", lwd=3, pch=7, axes=F, xlab='', ylab='', main='')
points(y, type='h')
axis(4)
print(names(x))
axis(1, at=1:length(x), labels=names(x))
par(op)
}
pareto(attenu$event)
title(main="Graphique de Paréto avec effectifs cumulés")
x <- apply(HairEyeColor, 2, sum) pie(x) title(main="Camembert")
On constate que les êtres humains perçoivent mal les différences entre des angles ou des surfaces. Or c'est justement ce sur quoi reposent les camemberts : ce n'est donc pas toujours un bon moyen de présenter l'information (c'est beau, mais peu informatif).
Quand on a beaucoup de valeurs, les "dotplots" peuvent remplacer les graphiques en colonnes.
x <- apply(HairEyeColor, 2, sum) dotchart(x)
Ils restent très lisibles.
data(Cars93) dotchart(table(Cars93$Manufacturer))
library(nlme) data(Milk) dotchart(table(Milk$Cow))
On a cette fois-ci deux séries de nombres, en parallèle : il s'agit le plus souvent d'un tableau à deux colonnes, chaque ligne du tableau (chaque couple de valeurs) correspondant à un individu.
data(cars)
plot(cars$dist ~ cars$speed,
xlab = "Speed (mph)",
ylab = "Stopping distance (ft)",
las = 1)
title(main = "Nuage de points")
On peut ajouter un diagramme de dispersion unidimensionnel dans les marges.
plot(cars$dist ~ cars$speed,
xlab = "Speed (mph)",
ylab = "Stopping distance (ft)",
las = 1)
title(main = "cars data")
rug(side=1, jitter(cars$speed, 5))
rug(side=2, jitter(cars$dist, 20))
On peut aussi ajouter une boite à moustache.
op <- par()
layout( matrix( c(2,1,0,3), 2, 2, byrow=T ),
c(1,6), c(4,1),
)
par(mar=c(1,1,5,2))
plot(cars$dist ~ cars$speed,
xlab='', ylab='',
las = 1)
rug(side=1, jitter(cars$speed, 5) )
rug(side=2, jitter(cars$dist, 20) )
title(main = "cars data")
par(mar=c(1,2,5,1))
boxplot(cars$dist, axes=F)
title(ylab='Stopping distance (ft)', line=0)
par(mar=c(5,1,1,2))
boxplot(cars$speed, horizontal=T, axes=F)
title(xlab='Speed (mph)', line=1)
par(op)
On peut tenter d'approcher ces données par une droite.
plot(dist ~ speed, data=cars) abline(lm(dist ~ speed, data=cars), col='red')
La fonction lowess permet d'avoir une courbe non nécessairement rectiligne (nous verrons d'autres exemples du même type quand nous parlerons de régression (linéaire, curvilinéaire, robuste)).
plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)",
las = 1)
lines(lowess(cars$speed, cars$dist, f = 2/3, iter = 3), col = "red")
title(main = "cars data")
Si l'une des variable est périodique, on peut la représenter sur un cercle.
On peut par exemple représenter l'affluence sur un site Web en fonction de l'heure du jour ou le volume des précipitation en fonction de la direction du vent. (La fonction stars permet aussi de faire ce genre de chose.)
x <- c(15, 9, 75, 90, 1, 1, 11, 5, 9, 8, 33, 11, 11,
20, 14, 13, 10, 28, 33, 21, 24, 25, 11, 33)
# j'ai essayé pendant près d'une demi-heure d'utiliser la commande stars
# pour faire ça, mais sans succès.
clock.plot <- function (x, col=rainbow(n), ...) {
if( min(x)<0 ) x <- x - min(x)
if( max(x)>1 ) x <- x/max(x)
n <- length(x)
if(is.null(names(x))) names(x) <- 0:(n-1)
m <- 1.05
plot(0, type='n', xlim=c(-m,m), ylim=c(-m,m),
axes=F, xlab='', ylab='', ...)
a <- pi/2 - 2*pi/200*0:200
polygon( cos(a), sin(a) )
v <- .02
a <- pi/2 - 2*pi/n*0:n
segments( (1+v)*cos(a), (1+v)*sin(a), (1-v)*cos(a), (1-v)*sin(a) )
segments( cos(a), sin(a), 0, 0, col='light grey', lty=3)
ca <- -2*pi/n*(0:50)/50
for (i in 1:n) {
a <- pi/2 - 2*pi/n*(i-1)
b <- pi/2 - 2*pi/n*i
polygon( c(0, x[i]*cos(a+ca), 0),
c(0, x[i]*sin(a+ca), 0),
col=col[i] )
v <- .1
text((1+v)*cos(a), (1+v)*sin(a), names(x)[i])
}
}
clock.plot(x, main="Affluence d'un site Web en fonction de l'heure")
On peut aussi découper l'ensemble des valeurs de Y en tranches et, pour chacune de ces tranches, faire une boite à moustaches, un histogramme ou une courbe de densité des valeurs de X correspondantes. Ces fonctions sont dans la bibliothèque lattice.
library(lattice) y <- cars$dist x <- cars$speed # vitesse <- shingle(x, co.intervals(x, number=6)) vitesse <- equal.count(x) histogram(~ y | vitesse)
bwplot(~ y | vitesse, layout=c(1,6))
densityplot(~ y | vitesse, aspect='xy')
densityplot(~ y | vitesse, layout=c(1,6))
Avant de connaître les graphiques en treillis, je découpais l'une des variables en quartiles, et je représentais la médiane (ou les quartiles) de l'autre pour chacun de ces quartiles.
y <- cars$dist
x <- cars$speed
q <- quantile(x)
o1 <- x<q[2]
o2 <- q[2]<x & x<q[3]
o3 <- q[3]<x & x<q[4]
o4 <- q[4]<x
dotchart(c(median(y[o1]), median(y[o2]),
median(y[o3]), median(y[o4])),
labels=as.character(1:4),
xlab="speed", ylab="distance")
my.dotchart <- function (y,x,...) {
x <- as.matrix(x)
z <- NULL
for (i in 1:dim(x)[2]) {
q <- quantile(x[,i])
for (j in 1:4) {
z <- append(z, median(y[ q[j] <= x[,i] & x[,i] <= q[j+1] ]))
names(z)[length(z)] <- paste(colnames(x)[i], as.character(j))
}
}
dotchart(z,...)
}
my.dotchart(y,x, xlab="speed", ylab="distance")
my.dotchart <- function (y,x,...) {
x <- as.matrix(x)
z <- NULL
for (i in 1:dim(x)[2]) {
q <- quantile(x[,i])
for (j in 1:4) {
ya <- y[ q[j] <= x[,i] & x[,i] <= q[j+1] ]
z <- rbind(z, quantile(ya))
rownames(z)[dim(z)[1]] <- paste(colnames(x)[i], as.character(j))
}
}
dotchart(t(z))
}
my.dotchart(y,x, xlab="speed", ylab="distance")
Je n'arrive pas à superposer ce genre de diagramme : faisons les choses à la main.
my.dotchart <- function (y,x,...) {
x <- as.matrix(x)
z <- NULL
for (i in 1:dim(x)[2]) {
q <- quantile(x[,i])
for (j in 1:4) {
ya <- y[ q[j] <= x[,i] & x[,i] <= q[j+1] ]
z <- rbind(z, quantile(ya))
rownames(z)[dim(z)[1]] <- paste(colnames(x)[i], as.character(j))
}
}
xmax <- max(z)
xmin <- min(z)
n <- dim(z)[1]
plot( z[,3], 1:n, xlim=c(xmin,xmax), ylim=c(1,n), axes=F, frame.plot=T, pch='.', ... )
axis(1)
axis(2, at=1:n, las=1)
abline( h=1:n, lty=3 )
# médiane
points( z[,3], 1:n, pch=16, cex=3 )
# quartiles
segments( z[,2], 1:n, z[,4], 1:n, lwd=7 )
# min et max
segments( z[,1], 1:n, z[,5], 1:n )
}
my.dotchart(y,x, xlab="speed", ylab="distance")
On vient en fait de réinventer la boite à moustaches...
On peut faire une boite à moustaches par valeur de la variable qualitative (ici, on compare l'efficacité de différents insecticides).
data(InsectSprays)
boxplot(count ~ spray, data = InsectSprays,
xlab = "Type of spray", ylab = "Insect count",
main = "InsectSprays data", varwidth = TRUE, col = "lightgray")
Les violons sont une variante des boites à moustaches, qui représentent aussi les densités. On peut ainsi repérer plus facilement des densités bimodales.
library(Simple) n <- 1000 k <- 10 x <- factor(sample(1:5, n, replace=T)) m <- rnorm(k,sd=2) s <- sample( c(rep(1,k-1),2) ) y <- rnorm(n, m[x], s[x]) x[x %in% order(m)[c(1,5)]] <- order(m)[1] simple.violinplot(y~x, col='pink')
On peut reprendre la fonction my.dotchart, définie plus haut, et la modifier pour qu'elle accepte aussi les variables qualitatives.
my.dotchart <- function (y,x,...) {
x <- data.frame(x)
z <- NULL
cn <- NULL
for (i in 1:dim(x)[2]) {
if( is.numeric(x[,i]) ) {
q <- quantile(x[,i])
for (j in 1:4) {
ya <- y[ q[j] <= x[,i] & x[,i] <= q[j+1] ]
z <- rbind(z, quantile(ya))
cn <- append(cn, paste(colnames(x)[i], as.character(j)))
}
} else {
for (j in levels(x[,i])) {
ya <- y[ x[,i] == j ]
z <- rbind(z, quantile(ya))
cn <- append(cn, paste(colnames(x)[i], as.character(j)))
}
}
}
xmax <- max(z)
xmin <- min(z)
n <- dim(z)[1]
plot( z[,3], 1:n, xlim=c(xmin,xmax), ylim=c(1,n), axes=F, frame.plot=T, pch='.', ... )
axis(1)
axis(2, at=1:n, labels=cn, las=1)
abline( h=1:n, lty=3 )
# médiane
points( z[,3], 1:n, pch=16, cex=3 )
# quartiles
segments( z[,2], 1:n, z[,4], 1:n, lwd=7 )
# min et max
segments( z[,1], 1:n, z[,5], 1:n )
}
spray <- InsectSprays$spray
y <- InsectSprays$count
my.dotchart(y,spray, xlab="count", ylab="spray")
stripchart(InsectSprays$count ~ InsectSprays$spray, method='jitter')
On aimerait faire un « stripchart » avec différentes couleurs pour distinguer les classes, un peu comme on sait faire en dimensions supérieures.
data(iris)
plot(iris[1:4], pch=21, bg=c("red", "green", "blue")[as.numeric(iris$Species)])
Le problème, c'est que la commande stripchart n'accepte pas cet argument « bg ». On fait donc la même chose à la main. En fait, ça n'était pas une bonne idée : on ne voit rien...
a <- InsectSprays$count
b <- rnorm(length(a))
plot(b~a, pch=21,
bg=c("red", "green", "blue", "cyan", "yellow", "black")
[as.numeric(InsectSprays$spray)])
a <- as.vector(t(iris[1]))
b <- rnorm(length(a))
plot(b~a, pch=21, bg=c("red", "green", "blue")[as.numeric(iris$Species)])
Regardons sur des données construites à cet effet.
do.it <- function (v) {
n <- 100
y <- sample( 1:3, n, replace=T )
a <- runif(1)
b <- runif(1)
c <- runif(1)
x <- ifelse( y==1, a+v*rnorm(n),
ifelse( y==2, b+v*rnorm(n), c+v*rnorm(n) ))
r <- rnorm(n)
plot( r ~ x, pch=21, bg=c('red', 'green', 'blue')[y] )
}
do.it(.1)
Il faut vraiment que les groupes soient bien séparés pour y voir quelque chose.
do.it(.05)
On peut aussi faire plusieurs histogrammes.
hists <- function (x, y, ...) {
y <- factor(y)
n <- length(levels(y))
op <- par( mfcol=c(n,1) )
b <- hist(x, ..., plot=F)$breaks
for (l in levels(y)){
hist(x[y==l], breaks=b, probability=T, ylim=c(0,.3), main=l, col='lightblue', ...)
points(density(x[y==l]), type='l', lwd=3, col='red')
}
par(op)
}
hists(InsectSprays$count, InsectSprays$spray)
On peut aussi faire ça avec des treillis.
library(lattice) histogram( ~ count | spray, data=InsectSprays)
densityplot( ~ count | spray, data = InsectSprays )
bwplot( ~ count | spray, data = InsectSprays )
bwplot( ~ count | spray, data = InsectSprays, layout=c(1,6) )
On peut calculer un « rapport de corrélation » : si X est la variable qualitative et Y la variable quantitative, on peut écrire la variance de Y comme somme d'une « variance expliquée par » et d'une « variance résiduelle » ; le rapport de corrélation est le quotient (parfois on prend la racine de ce quotient, parfois non...) de la variance expliquée et de la variance de Y.
> aov( count ~ spray, data = InsectSprays )
Call:
aov(formula = r)
Terms:
spray Residuals
Sum of Squares 2668.833 1015.167
Deg. of Freedom 5 66
Residual standard error: 3.921902
Estimated effects may be unbalancedIci, la variance expliquée est 2669, la variance résiduelle est 1015. Le rapport de corrélation est donc (quotient de la variance expliquée et de la variance totale),
> s <- summary(aov( count ~ spray, data = InsectSprays ))[[1]][,2] > s[1]/(s[1]+s[2]) [1] 0.7244
On résume cela en disant que le type d'insecticide explique 70% de la variance du nombre d'insectes.
On peut avoir cette valeur autrement (au milieu de plein d'autres choses, sur lesquelles nous reviendrons quand nous aborderons la régression et l'anova (analyse de la variance)).
> summary( lm( count ~ spray, data = InsectSprays ) )
Call:
lm(formula = count ~ spray, data = InsectSprays)
Residuals:
Min 1Q Median 3Q Max
-8.333 -1.958 -0.500 1.667 9.333
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.5000 1.1322 12.807 < 2e-16 ***
sprayB 0.8333 1.6011 0.520 0.604
sprayC -12.4167 1.6011 -7.755 7.27e-11 ***
sprayD -9.5833 1.6011 -5.985 9.82e-08 ***
sprayE -11.0000 1.6011 -6.870 2.75e-09 ***
sprayF 2.1667 1.6011 1.353 0.181
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 3.922 on 66 degrees of freedom
Multiple R-Squared: 0.7244, Adjusted R-squared: 0.7036
F-statistic: 34.7 on 5 and 66 DF, p-value: < 2.2e-16
On peut présenter ces données sous forme d'un tableau, les lignes représentant les différentes valeurs de la première variable, les colonnes celles de la seconde, les cases du tableau contiennent les effectifs.
data(HairEyeColor) a <- as.table( apply(HairEyeColor, c(1,2), sum) ) barplot(a, legend.text = attr(a, "dimnames")$Hair)
La table a contient :
Eye
Hair Brown Blue Hazel Green
Black 68 20 15 5
Brown 119 84 54 29
Red 26 17 14 14
Blond 7 94 10 16Voici d'autres manières de présenter la même chose.
barplot(a, beside=TRUE, legend.text = attr(a, "dimnames")$Hair)
barplot(t(a), legend.text = attr(a, "dimnames")$Eye)
barplot(t(a), beside=TRUE, legend.text = attr(a, "dimnames")$Eye)
On peut remplacer les lignes par les « profils-lignes », i.e., des fréquences marginales (on peut aussi faire ça avec les colonnes). On peut ensuite représenter chaque ligne par une barre.
> b <- a / apply(a, 1, sum)
> apply(b, 1, sum)
Black Brown Red Blond
1 1 1 1
> options(digits=2)
> b
Eye
Hair Brown Blue Hazel Green
Black 0.630 0.19 0.139 0.046
Brown 0.416 0.29 0.189 0.101
Red 0.366 0.24 0.197 0.197
Blond 0.055 0.74 0.079 0.126
> c <- t( t(a) / apply(a, 2, sum) )
> apply(c, 2, sum)
Brown Blue Hazel Green
1 1 1 1
> c
Eye
Hair Brown Blue Hazel Green
Black 0.309 0.093 0.16 0.078
Brown 0.541 0.391 0.58 0.453
Red 0.118 0.079 0.15 0.219
Blond 0.032 0.437 0.11 0.250
b <- a / apply(a, 1, sum)
barplot(t(b))
c <- t( t(a) / apply(a, 2, sum) ) barplot(c)
C'est déjà ce que fait la commande mosaicplot (en plus, la largeur des barres dépend des effectifs marginaux).
plot(a)
plot(t(a))
On peut ajouter des couleurs :
plot(a, col=heat.colors(dim(a)[2]))
ou plus simplement :
plot(a, color=T)
Je n'ai pas réussi à avoir simplement à la fois des couleurs (pour distinguer les valeurs de la première variable) et des hachures (pour la seconde).
On peut aussi demander à R d'indiquer quelles sont les classes dont les effectifs sont significativement importants ou faibles.
plot(a, shade=T)
plot(t(a), shade=T)
On peut représenter chaque ligne du tableau par une courbe (l'axe des abscisses est alors utilisé pour l'autre variable qualitative et l'axe vertical représente les effectifs).
data(HairEyeColor)
a <- apply(HairEyeColor, c(1,2) , sum)
qualplot <- function (a) {
matplot( row(a), a, type='l', axes=F,
col = 1:dim(a)[2]+1,
lty = 1:dim(a)[2],
lwd=3,
xlab=names(dimnames(a))[1], ylab=names(dimnames(a))[2] )
axis(1, 1:dim(a)[1], row.names(a))
axis(2)
legend(1, max(a), row.names(t(a)),
lwd=3, cex=1.5,
col=1:dim(a)[2]+1, lty=1:dim(a)[2])
}
# Pour une utilisation interactive
qualplots <- function (a) {
op = par(ask=T)
qualplot(a)
qualplot(t(a))
par(op)
}
qualplot(a)
qualplot(t(a))
On peut préalablement modifier le tableau pour que les effectifs marginaux soient égaux, i.e., remplacer le tableau par les profils-lignes ou par les profils-colonnes.
qualplotfreq <- function (a) {
a <- t( t(a) / apply(a,2,sum) )
qualplot(a)
}
qualplotsfreq <- function (a) {
op = par(ask=T)
qualplotfreq(a)
qualplotfreq(t(a))
par(op)
}
qualplotfreq(a)
qualplotfreq(t(a))
Si les variables n'ont que deux valeurs, on peut aussi utiliser la commande fourfoldplot.
data(bacteria) fourfoldplot( table(bacteria$y, bacteria$ap) )
La superficie des quarts de cercle est proportiionnelle aux effectifs. Les deux variables sont liées si les quarts de cercle opposés (ceux de même couleur) ont des tailles comparables et différentes de celles des deux autres quarts de cercle. Pour faciliter la comparaison, on trace aussi les intervalles de confiance à 95% : s'ils ne se chevauchent pas, la différence est statistiquement significative. Dans l'exemple précédent, la différence est significative. Dans les exemples suivants, les variables sont liées pour la semaine 6, mais pas pour les autres (ou du moins, pas de manière statistiquement significative).
fourfoldplot( table(bacteria$y, bacteria$ap, bacteria$week) )
On peut représenter trois variables quantitatives par des disques dans le plan : les deux premières variables fournissent les coordonnées du centre, la troisième le diamètre (si on veut montrer que les variations de cette variable sont fortes) ou l'aire (si on veut montrer qu'elles sont faibles, ou si on est honnète).
n <- 50
x <- rnorm(n)
y <- rnorm(n)
z <- rnorm(n)
my.renorm <- function (z) {
z <- abs(z)
z <- 10*z/max(z)
z
}
z <- my.renorm(z)
plot(x,y,cex=z)
plot(x,y,cex=z,pch=16,col='red') points(x,y,cex=z)
Il y a des variantes : on peut ajouter une variable (le plus souvent qualitative) en faisant varier la couleur des disques ;
u <- sample(c('red','green','blue'),n,replace=T)
plot(x,y,cex=z,pch=16,col=u)
points(x,y,cex=z)
on peut ajouter plusieurs variables quantitatives en remplaçant le disque par des cercles concentriques, etc.
z2 <- rnorm(n) z2 <- my.renorm(z2) plot(x,y,cex=z) points(x,y,cex=z2,col='red')
# Autre renormalisation (pertinente s'il n'y a pas de zéro)
my.renorm <- function (z) {
z <- (z-min(z)) / (max(z)-min(z))
z <- 1+9*z
z
}
z <- my.renorm(z)
z2 <- my.renorm(z2)
plot(x,y,cex=z)
points(x,y,cex=z2,col='red')
On peut aussi remplacer les cercles concentriques par des diagrammes étoilés.
n <- 50
x <- runif(n)
y <- runif(n)
z1 <- rnorm(n)
z2 <- rnorm(n)
z3 <- rnorm(n)
z4 <- rnorm(n)
z5 <- rnorm(n)
stars( data.frame(z1,z2,z3,z4,z5), location=cbind(x,y),
labels=NULL, len=1/sqrt(n)/2 )
v <- .2
n <- 50
x <- runif(n)
y <- runif(n)
z1 <- x+y+v*rnorm(n)
z2 <- x*y+v*rnorm(n)
z3 <- x^2 + y^2 + v*rlnorm(n)
stars( data.frame(z1,z2,z3), location=cbind(x,y),
labels=NULL, len=1/sqrt(n)/2, axes=T,
draw.segments=T, col.segments=1:5 )
Il s'agit de tracer plusieurs variables Y1, Y2, etc. en fonction de X. On pourrait superposer les courbes les unes aux autres, avec la commande matplot.
n <- 10
d <- data.frame(y1=abs(rnorm(n)),
y2=abs(rnorm(n)),
y3=abs(rnorm(n)),
y4=abs(rnorm(n)),
y5=abs(rnorm(n))
)
matplot(d, type='l')
On pourrait utiliser un diagramme en barres.
barplot(t(as.matrix(d)))
Pour faire un "line chart", on procède ainsi. On commence par tracer Y1 en fonction de X, on colorie en dessous de la courbe, on trace Y1+Y2 en fonction de X, on colorie entre les deux courbes, etc.
line.chart <- function (d) {
m <- d
m <- t(apply(m,1,cumsum))
print(m)
n1 <- dim(m)[1]
n2 <- dim(m)[2]
col = rainbow(n)
plot.new()
plot.window(xlim=c(1,n1), ylim=c(min(m),max(m)))
axis(1)
axis(2)
for (i in n2:1) {
polygon(c(1:n1,n1,1), c(m[,i],0,0), col=col[i], border=0)
}
for (i in n2:1) {
lines(m[,i], lwd=2)
}
}
line.chart(d)
On peut s'amuser à ajouter différents effets 3D -- mais R n'est pas un logiciel de dessin...
Un camembert avec une ombre.
shaded.pie <- function (...) {
pie(...)
op <- par(new=T)
a <- seq(0,2*pi,length=100)
for (i in (256:64)/256) {
r <- .8-.1*(1-i)
polygon( .1+r*cos(a), -.2+r*sin(a), border=NA, col=rgb(i,i,i))
}
par(new=T)
pie(...)
par(op)
}
x <- rpois(10,5)
x <- x[x>0]
shaded.pie(x)
Des barres en 3D (avec kchart : je laisse au lecteur le soin de faire pareil avec R).
On pourrait aussi mettre les barres les unes derrière les autres, on pourrait aussi faire des rubans en 3D.
On pourrait aussi déléguer la tache de dessin à un logiciel comme PoVRay...
povray.hist <- function (x, ...) {
x <- hist(x, plot=F, ...)$counts
if( min(x)<0 ) x <- x-min(x)
x <- x/max(x)*1.5
n <- length(x)
w <- 2/n
e <- .01
my.cat <- function (s) { cat(s); cat("\n") }
my.cat('#include "colors.inc"')
my.cat('#include "stones.inc"')
my.cat('#include "woods.inc"')
my.cat('#declare rd=seed(0);')
my.cat('#declare T1 = texture { T_Stone16 scale .2}')
my.cat('#declare T2 = texture { T_Wood9 }')
my.cat('background { color White }')
my.cat('camera{')
my.cat(' location <.6,2,-1.5>')
my.cat(' look_at <0,.6,0>')
my.cat('}')
my.cat('light_source{<-1, 4, -2> color .5*White}')
my.cat('light_source{<+1, 4, -2> color White}')
my.cat('light_source{<1,.5,-2> color White}')
my.cat('// Ground ')
my.cat('#declare e = .1 ;')
my.cat('box{ <-1-e, -.02, -.25-e>, <1+e, 0, .25+e> texture {T1} }')
my.cat('// Boxes')
for (i in 1:n) {
my.cat(paste( 'box{<', -1+(i-1)*w +e, ',', x[i], ',', -w/2, '>, <',
-1+i*w, ',', 0, ',', w/2, '> ',
'texture{T2 rotate rand(rd)*90*x ',
'translate 2*<rand(rd), rand(rd), rand(rd)> }}' ))
}
}
povray.hist(faithful$eruptions)
Une animation, c'est simplement une suite d'images statiques : l'exemple suivant représente le graphique quantiles-quantiles de x^a, pour différentes valeurs de a. Dans un coin, on a mis l'histogramme de x^a, une estimation de sa densité et
n <- 500
set.seed(7236)
x <- rnorm(n,10)
N <- 50
a <- seq(.002,10,length=N)
for (i in factor(1:N)) {
png(filename=sprintf("animation_1_%03d.png", i))
y <- x^a[i]
qqnorm(y, main=paste("a =", signif(a[i])))
qqline(y, col='red')
par(fig=c(.7,.9,.2,.4),mar=c(0,0,0,0),new=T)
hist(y, col='light blue', probability=T,axes=F,main='')
lines(density(y), col='red', lwd=3)
curve(dnorm(x,mean(y),sd=sd(y)),col='blue',lty=2,lwd=3, add=T)
dev.off()
}On convertit cette suite d'images en une animation, par exemple avec ImageMagick.
convert -delay 50 -loop 50 *.png animation_1.gif
Le format MNG est au PNG ce que le GIF animé est au GIF :
convert -delay 50 -loop 50 *.png animation_1.mng
A FAIRE
Pour avoir du divx ???
La commande pairs produit un tableau de nuages de points.
data(LifeCycleSavings) plot(LifeCycleSavings)
On peut aussi demander le tableau des corrélations.
> cor(LifeCycleSavings)
sr pop15 pop75 dpi ddpi
sr 1.00 -0.46 0.32 0.22 0.30
pop15 -0.46 1.00 -0.91 -0.76 -0.05
pop75 0.32 -0.91 1.00 0.79 0.03
dpi 0.22 -0.76 0.79 1.00 -0.13
ddpi 0.30 -0.05 0.03 -0.13 1.00On peut configurer le graphique, en mettant d'autres choses dans les cases : des histogrammes, des coefficients de corrélation, etc. (néanmoins, ce n'est pas suffisemment configurable à mon goût, car les fonctions que l'on passe en argument à la commande pairs ne connaissent pas le numéro de la ligne et de la colonne courrantes).
panel.hist <- function(x, ...) {
usr <- par("usr");
on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y, ...)
}
# Le coefficient de corrélation
my.panel.smooth <-
function (x, y, col = par("col"), bg = NA, pch = par("pch"),
cex = 1, col.smooth = "red", span = 2/3, iter = 3, ...)
{
points(x, y, pch = pch, col = col, bg = bg, cex = cex)
ok <- is.finite(x) & is.finite(y)
if (any(ok))
lines(lowess(x[ok], y[ok], f = span, iter = iter), col = col.smooth,
...)
usr <- par('usr')
text( (usr[1]+usr[2])/2, (usr[3]+9*usr[4])/10,
floor(100*cor(x,y))/100,
col='blue', cex=3, adj=c(.5,1) )
}
pairs(LifeCycleSavings,
diag.panel=panel.hist,
upper.panel=panel.smooth,
lower.panel=my.panel.smooth);
On peut aussi représenter les coefficients de corrélation par des couleurs. Dans l'exemple suivant, on voit que y est corrélée avec x1, x4 et x5 et que x4 et x5 sont corrélées. (Nous retrouverons ce genre d'exemple quand nous parlerons de régression multiple et de sélection des variables : face à de telles données, si on voulait prédire y en fonction des xi, on se limiterait aux prédictives x1 et x4 ou x1 et x5 (comme x4 et x5 sont fortement corrélées, une seule suffit)).
cor.plot <- function (x) {
n <- dim(x)[1]
m <- dim(x)[2]
N <- 1000
col = topo.colors(N)
plot(NA, xlim=c(0,1.2), ylim=c(-1,0))
for (i in 1:n) {
for (j in 1:m) {
polygon( c((j-1)/m, (j-1)/m, j/m, j/m),
-c((i-1)/m, i/m, i/m, (i-1)/m),
col = col[ N*(x[i,j]+1)/2 ] )
}
}
for (i in 1:N) {
polygon( c(1.1, 1.1, 1.2, 1.2),
-c((i-1)/N, i/N, i/N, (i-1)/N),
col = col[N-i+1],
border=NA )
}
# Exercice : ajouter une légende et le nom des variables
}
n <- 200
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)
x <- rnorm(n)
x4 <- x + .1*rnorm(n)
x5 <- x + .1*rnorm(n)
y <- 1 + x1 + x4 + rnorm(n)
d <- data.frame(y,x1,x2,x3,x4,x5)
cor.plot(cor(d))
On peut aussi faire une boite à moustache des valeurs de Y sur chaque quantile de chaque variable Xi (ça n'est pas symétrique : une variable joue une rôle particulier).
my.dotchart(LifeCycleSavings[,1], LifeCycleSavings[,-1],
xlab='savings', ylab='')
Avec des vraies boites à moustaches :
to.factor.vector <- function (x, number=4) {
resultat <- NULL
intervalles <- co.intervals(x,number,overlap=0)
for (i in 1:number) {
if( i==1 ) intervalles[i,1] = min(x)
else
intervalles[i,1] <- intervalles[i-1,2]
if( i==number )
intervalles[i,2] <- max(x)
}
for (valeur in x) {
r <- NA
for (i in 1:number) {
if( valeur >= intervalles[i,1] & valeur <= intervalles[i,2] )
r <- i
}
resultat <- append(resultat, r)
}
factor(resultat, levels=1:number)
}
to.factor <- function (x, number=4) {
if(is.vector(x))
r <- to.factor.vector(x, number)
else {
r <- NULL
for (v in x) {
a <- to.factor.vector(v)
if( is.null(r) )
r <- data.frame(a)
else
r <- data.frame(r,a)
}
names(r) <- names(x)
}
r
}
x <- to.factor(LifeCycleSavings[,-1])
y <- LifeCycleSavings[,1]
y <- as.vector(matrix(y,nr=length(y),ncol=dim(x)[2]))
for (i in names(x)) {
levels(x[[i]]) <- paste(i, levels(x[[i]]))
}
col <- gl( dim(x)[2], length(levels(x[,1])),
labels=rainbow( dim(x)[2] ))
col <- as.vector(col)
x <- factor(as.vector(as.matrix(x)))
boxplot(y~x, horizontal=T, las=1, col=col)
On doit aussi pouvoir faire ça avec des treillis (je n'ai pas réussi à mettre ça dans une boucle : rien ne s'affichait).
bwplot( ~ LifeCycleSavings[,1] | equal.count(LifeCycleSavings[,2], number=4),
layout=c(1,4) )
bwplot( ~ LifeCycleSavings[,1] | equal.count(LifeCycleSavings[,3], number=4),
layout=c(1,4) )
bwplot( ~ LifeCycleSavings[,1] | equal.count(LifeCycleSavings[,4], number=4),
layout=c(1,4) )
bwplot( ~ LifeCycleSavings[,1] | equal.count(LifeCycleSavings[,5], number=4),
layout=c(1,4) )
Si on dispose de beaucoup de variables (une demi-douzaine, quantitatives) et peu d'individus (une dizaine), on peut faire des graphiques étoilés.
data(mtcars)
stars(mtcars[, 1:7], key.loc = c(14, 2),
main = "Motor Trend Cars : stars(*, full = F)", full = FALSE)
stars(mtcars[, 1:7], key.loc = c(14, 1.5),
main = "Motor Trend Cars : full stars()",flip.labels=FALSE)
palette(rainbow(12, s = 0.6, v = 0.75))
stars(mtcars[, 1:7], len = 0.8, key.loc = c(12, 1.5),
main = "Motor Trend Cars", draw.segments = TRUE)
stars(mtcars[, 1:7], locations = c(0,0), radius = FALSE,
key.loc=c(0,0), main="Motor Trend Cars", lty = 2)
R sait tracer des surfaces définie par une équation de la forme z = f(x,y). Elles sont décrites par deux vecteurs -- les valeurs des deux premières coordonnées -- et un tableau -- les valeurs de la dernière coordonnée.
On peut représenter une telle surface sous forme de fil de fer
# Tiré du manuel
x <- seq(-10, 10, length=50)
y <- x
f <- function(x,y) {
r <- sqrt(x^2+y^2)
10 * sin(r)/r
}
z <- outer(x, y, f)
z[is.na(z)] <- 1
persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
shade=.5,
xlab = "X", ylab = "Y", zlab = "Z")
# Exemple du manuel data(volcano) z <- 2 * volcano # Exaggerate the relief x <- 10 * (1:nrow(z)) # 10 meter spacing (S to N) y <- 10 * (1:ncol(z)) # 10 meter spacing (E to W) persp(x, y, z, theta = 120, phi = 15, scale = FALSE, axes = FALSE) # Voir aussi les autres exemples dans # demo(persp)
de lignes de niveau
# Toujours un exemple du manuel
data("volcano")
rx <- range(x <- 10*1:nrow(volcano))
ry <- range(y <- 10*1:ncol(volcano))
ry <- ry + c(-1,1) * (diff(rx) - diff(ry))/2
tcol <- terrain.colors(12)
opar <- par(pty = "s", bg = "lightcyan")
plot(x = 0, y = 0,type = "n", xlim = rx, ylim = ry, xlab = "", ylab = "")
u <- par("usr")
rect(u[1], u[3], u[2], u[4], col = tcol[8], border = "red")
contour(x, y, volcano, col = tcol[2], lty = "solid", add = TRUE,
vfont = c("sans serif", "plain"))
title("A Topographic Map of Maunga Whau", font = 4)
abline(h = 200*0:4, v = 200*0:4, col = "lightgray", lty = 2, lwd = 0.1)
par(opar)
avec éventuellement des couleurs.
# Exemple du manuel data(volcano) x <- 10*(1:nrow(volcano)) y <- 10*(1:ncol(volcano)) image(x, y, volcano, col = terrain.colors(100), axes = FALSE) contour(x, y, volcano, levels = seq(90, 200, by=5), add = TRUE, col = "peru") axis(1, at = seq(100, 800, by = 100)) axis(2, at = seq(100, 600, by = 100)) box() title(main = "Maunga Whau Volcano", font.main = 4)
La bibliothèque lattice contient des fonctions wireframe et cloud.
data(volcano) x <- 10*(1:nrow(volcano)) x <- rep(x, ncol(volcano)) y <- 10*(1:ncol(volcano)) y <- rep(y, each=nrow(volcano)) z <- as.vector(volcano) wireframe( z ~ x * y )
cloud( z ~ x * y )
Pour ceux qui disposent des lunettes adéquates, on peut aussi faire des dessins stéréoscopiques.
data(iris)
print(cloud(Sepal.Length ~ Petal.Length * Petal.Width,
data = iris, cex = .8,
groups = Species,
subpanel = panel.superpose,
main = "Stereo",
screen = list(z = 20, x = -70, y = 3)),
split = c(1,1,2,1), more = TRUE)
print(cloud(Sepal.Length ~ Petal.Length * Petal.Width,
data = iris, cex = .8,
groups = Species,
subpanel = panel.superpose,
main = "Stereo",
screen = list(z = 20, x = -70, y = 0)),
split = c(2,1,2,1))
R lui-même ne sait pas faire d'animations, mais on peut exporter les données pour les observer dans un autre logiciel, par exemple xgobi. (Quand la fenêtre d'xgobi apparait, appuyez sur 'g' pour lancer l'animation.)
library(xgobi)
n <- 50
x <- seq(-10, 10, length=n)
y <- x
f <- function(x,y) {
r <- sqrt(x^2+y^2)
10 * sin(r)/r
}
z <- outer(x, y, f)
z[is.na(z)] <- 1
x <- matrix(x,nr=n,nc=n)
y <- matrix(y,nr=n,nc=n,byrow=T)
sombrero <- data.frame(x=as.vector(x),y=as.vector(y),z=as.vector(z))
xgobi(sombrero)
C'est aussi très utile pour observer des données statistiques dans des espaces de dimension plus grande.
http://www.public.iastate.edu/~dicook/JSS/paper/paper.html A FAIRE : on peut donner différents fichiers à xgobi : x.dat données x.row row names x.col col names x.glyphs x.colors mais quand on l'appelle depuis R ???
Quelques tutoriels :
http://industry.ebi.ac.uk/%257Ealan/VisWorkshop99/XGobi_Tutorial/index.html http://www.public.iastate.edu/~dicook/stat503.html
La nouvelle version de xgobi s'appelle ggobi :
http://www.ggobi.org/
Si les données sont assez complexes, on peut demander à xgobi de les afficher en même temps dans plusieurs fenêtres, sélectionner certains points dans l'une d'entre elle et voir où ils sont dans les autres (en anglais, ça s'appelle "brushing").
n <- 1000 x <- runif(3) + rnorm(n) y <- runif(3) + rnorm(n) z <- runif(3) + rnorm(n) t <- c(-5,-5,5) + rnorm(n) u <- c(-5,5,5) + rnorm(n) three.clusters <- data.frame(x,y,z,t,u) xgobi(three.clusters)
Xgobi permet aussi de faire des "graphes parallèles" : à chaque point correspond une ligne brisée, dont les ordonnées successives sont les différentes coordonnées du point. Très souvent, on ne voit rien dessus. Mais on peut demander à xgobi de ne tracer que certains des points et le graphe peut devenir lisible.
n <- 100 m <- matrix( rnorm(5*n)+c(1,-1,3,0,2), nr=n, nc=5, byrow=T ) matplot(1:5, t(m), type='l') title(main="Nuage de points homogène")
n <- 50 k <- 2 m <- matrix( rnorm(5*k*n)+runif(5*k,min=-10,max=10), nr=n, nc=5, byrow=T ) matplot(1:5, t(m), type='l') title(main="Nuage de points avec des (2) grumeaux")
n <- 50 k <- 5 m <- matrix( rnorm(5*k*n)+runif(5*k,min=-10,max=10), nr=n, nc=5, byrow=T ) matplot(1:5, t(m), type='l') title(main="Nuage de points avec des (5) grumeaux")
matplot(1:5, t(princomp(m)$scores), type='l') title(main="idem, après ACP")
matplot(1:5, t(m), type='l') title(main="Nuage de points avec des grumeaux moins visibles")
Il y a aussi une fonction parallel dans la bibliothèque lattice.
library(lattice) ?parallel
A FAIRE :
http://wsopuppenkiste.wiso.uni-goettingen.de/~dadler/rgl/
On dispose de trois variables quantitatives (très souvent, des pourcentages, en chimie).
library(MASS)
data(Skye)
ternary <- function(X, pch = par("pch"), lcex = 1,
add = FALSE, ord = 1:3, ...)
{
X <- as.matrix(X)
if(any(X) < 0) stop("X must be non-negative")
s <- drop(X %*% rep(1, ncol(X)))
if(any(s<=0)) stop("each row of X must have a positive sum")
if(max(abs(s-1)) > 1e-6) {
warning("row(s) of X will be rescaled")
X <- X / s
}
X <- X[, ord]
s3 <- sqrt(1/3)
if(!add)
{
oldpty <- par("pty")
on.exit(par(pty=oldpty))
par(pty="s")
plot(c(-s3, s3), c(0.5-s3, 0.5+s3), type="n", axes=FALSE,
xlab="", ylab="")
polygon(c(0, -s3, s3), c(1, 0, 0), density=0)
lab <- NULL
if(!is.null(dn <- dimnames(X))) lab <- dn[[2]]
if(length(lab) < 3) lab <- as.character(1:3)
eps <- 0.05 * lcex
text(c(0, s3+eps*0.7, -s3-eps*0.7),
c(1+eps, -0.1*eps, -0.1*eps), lab, cex=lcex)
}
points((X[,2] - X[,3])*s3, X[,1], ...)
}
ternary(Skye/100, ord=c(1,3,2))
Ici, les lignes du tableau ont une somme égale à 100.
> sum( apply(Skye,1,sum) != 100 ) [1] 0
Voici une autre manière de les tracer, d'après http://finzi.psych.upenn.edu/R/Rhelp01/archive/1000.html
tri <-
function(a, f, m, symb = 2, grid = F, ...)
{
ta <- paste(substitute(a))
tf <- paste(substitute(f))
tm <- paste(substitute(m))
tot <- 100/(a + f +m)
b <- f * tot
y <- b * .878
x <- m * tot + b/2
par(pty = "s")
oldcol <- par("col")
plot(x, y, axes = F, xlab = "", ylab = "",
xlim = c(-10, 110), ylim= c(-10, 110), type = "n", ...)
points(x,y,pch=symb)
par(col = oldcol)
trigrid(grid)
text(-5, -5, ta)
text(105, -5, tm)
text(50, 93, tf)
par(pty = "m")
invisible()
}
trigrid <-
function(grid = F)
{
lines(c(0, 50, 100, 0), c(0, 87.8, 0, 0)) #draw frame
if(!grid) {
for(i in 1:4 * 20) {
lines(c(i, i - 1), c(0, 2 * .878)) #side a-c (base)
lines(c(i, i + 1), c(0, 2 * .878))
T.j <- i/2 #side a-b (left)
lines(c(T.j, T.j + 2), c(i * .878, i * .878))
lines(c(T.j, T.j + 1), c(i * .878, (i - 2) * .878))
T.j <- 100 - i/2 #side b-c (right)
lines(c(T.j, T.j - 2), c(i * .878, i * .878))
lines(c(T.j, T.j - 1), c(i * .878, (i - 2) * .878))
}
} else {
for(i in 1:4 * 20) {
# draw dotted grid
lines(c(i, i/2), c(0, i * .878), lty = 4, col = 3)
lines(c(i, (50 + i/2)), c(0, .878 * (100 - i)), lty = 4, col = 3)
lines(c(i/2, (100 - i/2)), c(i * .878, i * .878), lty = 4, col = 3)
}
par(lty = 1, col = 1)
}
}
# some random data in three variables
c1<- runif(5, 10, 20)
c2<- runif(5, 1, 5)
c3 <- runif(5, 15, 25)
# basic plot
tri(c1,c2,c3)
# plot with different symbols and a grid tri(c1,c2,c3, symb=7, grid=T)
C'est un bon exercice d'écrire soi-même sa propre fonction. Par exemple, en plus de la grille, seriez-vous capable de remplir le triangle avec un dégradé (à chaque sommet on associe une couleur) ?
On peut ajouter une dimension, et utiliser un logiciel tel xgobi pour visualiser le résultat.
library(xgobi) data(PaulKAI) quadplot(PaulKAI, normalize = TRUE)
Il est aussi possible de transformer les données en son (c'est ce que font les dauphins pour se répérer, c'est ainsi que fonctionnent les radars, etc.)...
http://www.matthiasheymann.de/Download/Sonification.pdf http://www.matthiasheymann.de/Download/sound-manual.pdf
Interlude musical : certains s'amusent à transformer l'ADN en musique...
http://www.bbc.co.uk/radio4/science/thematerialworld_20031120.shtml http://www.dnamusiccentral.com/#About%20DNA
On peut représenter les différentes valeurs de la variable qualitative par des symboles ou des couleurs différentes.
data(iris)
plot(iris[1:4], pch=21,
bg=c("red", "green", "blue")[as.numeric(iris$Species)])
Si l'une des variables qualitative est binaire, on peut aussi utiliser la commande fourfoldplot, comme vu plus haut.
Quand on a juste une variable qualitative et une variable quantitative, on a vu qu'on pouvait faire des graphes à moustaches.
Avec deux variables qualitatives et une variable quantitative, ce n'est plus lisible, car les graphes à moustaches se superposent mal.
a <- rnorm(10)
b <- 1+ rnorm(10)
c <- 1+ rnorm(10)
d <- rnorm(10)
x <- c(a,b,c,d)
y <- factor(c( rep("A",20), rep("B",20)))
z <- factor(c( rep("U",10), rep("V",20), rep("U",10) ))
op = par(mfrow=c(2,2))
plot(x~y)
plot(x~z)
plot(x[z=="U"] ~ y[z=="U"], border="red", ylim=c(min(x),max(x)))
plot(x[z=="V"] ~ y[z=="V"], border="blue", add=T)
plot(x[y=="A"] ~ z[y=="A"], border="red", ylim=c(min(x),max(x)))
plot(x[y=="B"] ~ z[y=="B"], border="blue", add=T)
par(op)
On peut juste les mettres côte à côte (si les variables n'ont que deux valeurs, comme ici, ça reste lisible).
l <- rep("",length(x))
for (i in 1:length(x)){
l[i] <- paste(y[i],z[i])
}
l <- factor(l)
boxplot(x~l)
Par contre, au lieu de superposer des graphes à moustaches, on peut superposer des courbes correspondant à la moyenne des échantillons.
# l est une liste à deux éléments
myplot1 <- function (x, l, ...) {
t <- tapply(x,l,mean)
l1 <- levels(l[[1]])
l2 <- levels(l[[2]])
matplot(t,
type='l', lty=1, col=1:length(l2),
axes=F, ...)
axis(1, 1:2, l1)
axis(2)
lim <- par("usr")
legend(lim[1] + .05*(lim[2]-lim[1]), lim[4],
l2, lwd=1, lty=1, col=1:length(l2) )
}
op <- par(mfrow=c(1,2))
myplot1( x, list(y,z), ylim=c(0,2) )
myplot1( x, list(z,y), ylim=c(0,2) )
par(op)
Si le dessin n'est pas trop chargé, on peut tenter de l'enrichir en faisant figurer, en pointillés, les quartiles et les valeurs extrèmes de chaque échantillon.
# l est une liste à deux éléments
myplot3 <- function (x, l, ...) {
l1 <- levels(l[[1]])
l2 <- levels(l[[2]])
t0 <- tapply(x,l,min)
t1 <- tapply(x,l,function(x)quantile(x,.25))
t2 <- tapply(x,l,median)
t3 <- tapply(x,l,function(x)quantile(x,.75))
t4 <- tapply(x,l,max)
matplot(cbind(t0,t1,t2,t3,t4),
type='l',
lty=c(rep(3,length(l2)), rep(2,length(l2)),
rep(1,length(l2)), rep(2,length(l2)),
rep(3,length(l2)) ),
col=1:length(l2),
axes=F, ...)
axis(1, 1:2, l1)
axis(2)
lim <- par("usr")
legend(lim[1] + .05*(lim[2]-lim[1]), lim[4],
l2, lwd=1, lty=1, col=1:length(l2) )
}
op <- par(mfrow=c(1,2))
myplot3( x, list(y,z) )
myplot3( x, list(z,y) )
par(op)
http://cm.bell-labs.com/cm/ms/departments/sia/project/trellis/
Les graphiques en treillis consistent à dessiner un nuage de points en 3 dimensions (ou 4, ou 5), à découper l'espace en tranches, et à projeter ces tranches sur des sous-espaces de dimension 2.
On peut aussi les voirs comme des graphiques dans lesquels on sépare les individus en plusieurs groupes (les tranches ci-dessus ou les différentes valeurs d'une (ou plusieurs) variable(s) qualitative(s)) pour faire un graphique pour chacun d'entre eux. On se retrouve alors aves plusieurs graphiques, et ne pas un seul.
library(help=lattice) library(lattice) ?Lattice ?xyplot
L'exemple suivant (trois variables quantitatives) représente les épicentres de tremblements de terre : on les a représentés dans un espace de dimension 3, que l'on coupe en tranches.
data(quakes) Depth <- equal.count(quakes$depth, number=8, overlap=.1) xyplot(lat ~ long | Depth, data = quakes)
C'est beaucoup plus lisible que la projection sur l'un des plans de coordonnées (car on oublie alors l'une des variables), que les trois projections, ou que la projection sur un plan "judicieusement choisi" (nous verrons plus loin comment il est choisi : c'est l'analyse en composantes principales).
plot(lat ~ long, data=quakes)
op <- par(mfrow=c(2,2)) plot(lat ~ long, data=quakes) plot(lat ~ -depth, data=quakes) plot(-depth ~ long, data=quakes) par(op)
library(mva) biplot(princomp(quakes[1:3]))
pairs( princomp(quakes[1:3])$scores )
Un exemple avec une variable quantitative et 3 variables qualitatives, dans lequel, si on fixe les variables quantitatives, il n'y a qu'un seul individu.
data(barley) barchart(yield ~ variety | year * site, data=barley)
L'argument scales permet de modifier l'apparence des axes et de leurs graduations (ici, on évite que les lables se chevauchent).
barchart(yield ~ variety | year * site, data = barley,
ylab = "Barley Yield (bushels/acre)",
scales = list(x = list(0, abbreviate = TRUE,
minlength = 5)))
# Raté... dotplot(yield ~ variety | year * site, data = barley)
dotplot(variety ~ yield | year * site, data = barley)
On représente les différentes valeurs d'une des variables qualitatives par des couleurs ou des symboles différents (et on remarque que pour l'une des fermes, le rendement a augmenté au lieu de diminuer : en fait, les données des deux années ont été interverties), à l'aide de l'argument groups.
data(barley)
dotplot(variety ~ yield | site, groups = year,
data = barley,
layout = c(1, 6), aspect = .5, pch = 16,
col.line = c("grey", "transparent"),
panel = "panel.superpose",
panel.groups = "panel.dotplot")
Un exemple avec une variable quantitative, des variables qualitatives, mais dans lequel il y a plusieurs individus dans chaque classe.
library(nlme) data(bdf) d <- data.frame( iq=bdf$IQ.perf, sex=bdf$sex, den=bdf$denomina ) d <- d[1:100,] bwplot( ~ d$iq | d$sex + d$den )
histogram( ~ d$iq | d$sex + d$den )
densityplot( ~ d$iq | d$sex + d$den )
stripplot( ~ d$iq | d$sex + d$den )
Autre exemple : une variable quantitative en fonction d'une autre variable quantitative et d'une variable qualitative. (On transforme la variable quantitative en shingle à l'aide de la commande equal.count.)
d <- data.frame( x=bdf$aritPOST, y=bdf$sex, z=equal.count(bdf$langPRET) ) bwplot( ~ x | y + z, data=d )
histogram( ~ x | y + z, data=d )
densityplot( ~ x | y + z, data=d )
Un exemple avec uniquement des variables qualitatives (on calcule préalablement les effectifs, à l'aide de la commande table).
d <- data.frame( x= (bdf$IQ.perf>11), y=bdf$sex, z=bdf$denomina ) d <- as.data.frame(table(d)) barchart( Freq ~ x | y * z, data=d )
Les formules utilisées par la bibliothèque lattice ont, de prime abord, un air un peu compliqué.
x ~ y Représente x en fonction de y (un seul dessin)
x ~ y | z Représente x en fonction de y après avoir découpé
les données selon les valeurs de z
x ~ y | z1 * z2 Idem, selon les valeurs de (z1,z2)
x~y|z, groups=t Idem, mais on utilise un symbole différent (ou une
couleur différente) selon les valeurs de tLes exemples suivants représentent des données univariées dans chacune des cases (sauf le dernier).
~ y ~ y | z ~ y | z1 * z2 ~ y | z1 = z2, groups=t
Si la variable utilisée pour découper l'ensemble des observations en morceaux est qualitative, le découpage se fait selon ces valeurs, mais si elle est quantitative, il faut préalablement la convertir en une variable qualitative à l'aide (par exemple) de la fonction equal.count.
A un niveau très élémentaire, la bibliothèque lattice permet de positionner précisément plusieurs graphiques.
n <- 200 x <- rnorm(n) y <- x^3+rnorm(n) plot1 <- xyplot(y~x) plot2 <- bwplot(x) # Attention : l'ordre est xmin, ymin, xmax, ymax print(plot1, position=c(0,.2,1,1), more=T) print(plot2, position=c(0,0,1,.2), more=F)
Mais on pouvait déjà faire ça avec par(fig=...)
n <- 200 x <- rnorm(n) y <- x^4+rnorm(n) k <- .7 op <- par(mar=c(0,0,0,0)) # Attention : l'ordre est xmin, xmax, ymin, ymax par(fig=c(0,k,0,k)) plot(y~x) par(fig=c(0,k,k,1), new=T) boxplot(x, horizontal=T) par(fig=c(k,1,0,k), new=T) boxplot(y, horizontal=F) par(op)
A FAIRE : Mettre cet exemple un peu plus haut, quand je parlais de par() ?
En fait, les commandes xyplot, dotplot, histogram, densityplot, stripplot se contentent de découper les données en tranches et délèguent le tracé de chacun des sous-graphiques à d'autres fonctions (panel.xyplot, etc.), que l'on peut modifier. Dans ces fonctions, il ne faut pas utiliser les commandes usuelles de tracé (points, lignes, segments, polygon, abline, etc.), mais celles de la bibliothèque grid (lpoints, llines, lsegments, , panel.abline, etc.)
Voici l'un des exemples du manuel.
dotplot(variety ~ yield | site, data = barley, groups = year,
panel = function(x, y, subscripts, ...) {
dot.line <- trellis.par.get("dot.line")
panel.abline(h = y, col = dot.line$col,
lty = dot.line$lty)
panel.superpose(x, y, subscripts, ...)
},
key = list(space="right", transparent = TRUE,
points=list(pch=trellis.par.get("superpose.symbol")$pch[1:2],
col=trellis.par.get("superpose.symbol")$col[1:2]),
text=list(c("1932", "1931"))),
xlab = "Barley Yield (bushels/acre) ",
aspect=0.5, layout = c(1,6), ylab=NULL)
La fonction treillis.par.get est l'équivalent de la fonction par pour les graphiques standard : elle permet d'accéder aux paramètres graphiques. On peut les modifier à l'aide de la commande treillis.par.set.
> str(trellis.par.get()) List of 23 $ background :List of 1 ..$ col: chr "#909090" $ add.line :List of 3 ..$ col: chr "#000000" ..$ lty: num 1 ..$ lwd: num 1 $ add.text :List of 3 ..$ cex : num 1 ..$ col : chr "#000000" ..$ font: num 1 $ bar.fill :List of 1 ..$ col: chr "#00ffff" $ box.dot :List of 4 ..$ col : chr "#000000" ..$ cex : num 1 ..$ font: num 1 ..$ pch : num 16 $ box.rectangle :List of 3 ..$ col: chr "#00ffff" ..$ lty: num 1 ..$ lwd: num 1 $ box.umbrella :List of 3 ..$ col: chr "#00ffff" ..$ lty: num 2 ..$ lwd: num 1 $ dot.line :List of 3 ..$ col: chr "#aaaaaa" ..$ lty: num 1 ..$ lwd: num 1 $ dot.symbol :List of 4 ..$ cex : num 0.8 ..$ col : chr "#00ffff" ..$ font: num 1 ..$ pch : num 16 $ plot.line :List of 3 ..$ col: chr "#00ffff" ..$ lty: num 1 ..$ lwd: num 1 $ plot.symbol :List of 4 ..$ cex : num 0.8 ..$ col : chr "#00ffff" ..$ font: num 1 ..$ pch : num 1 $ reference.line :List of 3 ..$ col: chr "#aaaaaa" ..$ lty: num 1 ..$ lwd: num 1 $ strip.background:List of 1 ..$ col: chr [1:7] "#ffd18f" "#c8ffc8" "#c6ffff" "#a9e2ff" ... $ strip.shingle :List of 1 ..$ col: chr [1:7] "#ff7f00" "#00ff00" "#00ffff" "#007eff" ... $ superpose.line :List of 3 ..$ col: chr [1:7] "#00ffff" "#ff00ff" "#00ff00" "#ff7f00" ... ..$ lty: num [1:7] 1 1 1 1 1 1 1 ..$ lwd: num [1:7] 1 1 1 1 1 1 1 $ regions :List of 1 ..$ col: chr [1:100] "#FF80FF" "#FF82FF" "#FF85FF" "#FF87FF" ... $ superpose.symbol:List of 4 ..$ cex : num [1:7] 0.8 0.8 0.8 0.8 0.8 0.8 0.8 ..$ col : chr [1:7] "#00ffff" "#ff00ff" "#00ff00" "#ff7f00" ... ..$ font: num [1:7] 1 1 1 1 1 1 1 ..$ pch : chr [1:7] "o" "o" "o" "o" ... $ axis.line :List of 4 ..$ line: num 0 ..$ col : chr "#000000" ..$ lty : num 1 ..$ lwd : num 1 $ box.3d :List of 3 ..$ col: chr "#000000" ..$ lty: num 1 ..$ lwd: num 1 $ par.xlab.text :List of 3 ..$ cex : num 1 ..$ col : chr "#000000" ..$ font: num 1 $ par.ylab.text :List of 3 ..$ cex : num 1 ..$ col : chr "#000000" ..$ font: num 1 $ par.main.text :List of 3 ..$ cex : num 1.2 ..$ col : chr "#000000" ..$ font: num 2 $ par.sub.text :List of 3 ..$ cex : num 1 ..$ col : chr "#000000" ..$ font: num 2
La commande show.setings permet de voir à quoi tout cela correspond.
show.settings()
La fonction panel.abline trace des droites (ici, on donne un argument h, donc il s'agit de lignes horizontales, mais on pourrait aussi tracer des lignes verticales (argument v), des droites données par leur équation ou par le résultat d'une régression.
La fonction panel.superpose permet d'utiliser différentes couleurs ou différents symboles selon les valeur d'une des variables.
L'argument key définit la légende.
Voici un autre exemple : on superpose histogramme, densité et densité normale.
y <- cars$dist
x <- cars$speed
vitesse <- shingle(x, co.intervals(x, number=6))
histogram(~ x | vitesse, type = "density",
panel = function(x, ...) {
ps <- trellis.par.get('plot.symbol')
nps <- ps
nps$cex <- 1
trellis.par.set('plot.symbol', nps)
panel.histogram(x, ...)
panel.densityplot(x, col = 'brown', lwd=3)
panel.xyplot(x = jitter(x), y = rep(0, length(x)), col='brown' )
panel.mathdensity(dmath = dnorm,
args = list(mean=mean(x),sd=sd(x)),
lwd=3, lty=2, col='white')
trellis.par.set('plot.symbol', ps)
})
Voici une liste des fonctions (de haut niveau) qu'on peut avoir envie d'utiliser.
> apropos("^panel\..*")
[1] "panel.3dscatter" "panel.3dscatter.new" "panel.3dwire"
[4] "panel.abline" "panel.barchart" "panel.bwplot"
[7] "panel.cloud" "panel.densityplot" "panel.dotplot"
[10] "panel.fill" "panel.grid" "panel.histogram"
[13] "panel.levelplot" "panel.linejoin" "panel.lmline"
[16] "panel.loess" "panel.mathdensity" "panel.pairs"
[19] "panel.parallel" "panel.qq" "panel.qqmath"
[22] "panel.qqmathline" "panel.splom" "panel.stripplot"
[25] "panel.superpose" "panel.superpose.2" "panel.tmd"
[28] "panel.wireframe" "panel.xyplot" "panel.smooth"
> apropos("^prepanel\..*")
[1] "prepanel.default.bwplot" "prepanel.default.cloud"
[3] "prepanel.default.densityplot" "prepanel.default.histogram"
[5] "prepanel.default.levelplot" "prepanel.default.parallel"
[7] "prepanel.default.qq" "prepanel.default.qqmath"
[9] "prepanel.default.splom" "prepanel.default.tmd"
[11] "prepanel.default.xyplot" "prepanel.lmline"
[13] "prepanel.loess" "prepanel.qqmathline"La commande panel.superpose est un remplacement de panel.xyplot s'il y a un argument "groups". Je crois que les autres sont directement compréhensibles.
panel.smooth trace les points (comme panel.xyplot) et ajoute
une courbe de régression
Quand on trace des courbes, on constate qu'on parvient mieux à comparer la pente des tangentes (i.e., à dire dans quel graphique la variable augmente le plus vite) si cette pente est voisine de 45 degrés. On qualifie de "banking" le choix des unités permettant de se rapprocher de cette pente.
C'est particulièrement utile pour les séries temporelles : dans la deuxième figure, il saute aux yeux que la décroissance est plus lente. On peut aussi le voir dans la première figure, mais il faut regarder attentivement.
## banking
data(sunspot)
xyplot(sunspot ~ 1:37 ,type = "l",
scales = list(y = list(log = TRUE)),
sub = "log scales")
xyplot(sunspot ~ 1:37 ,type = "l", aspect="xy",
scales = list(y = list(log = TRUE)),
sub = "log scales")
Exercice : prendre le code d'une de ces fonctions, par exemple xyplot ou splom et le commenter. (Attention, c'est long!)
A FAIRE
Vincent Zoonekynd
<zoonek@math.jussieu.fr>
latest modification on Wed Oct 13 22:32:38 BST 2004