Histogramme et estimateur de densité

L'histogramme est un outil très utile pour visualiser la répartition des valeurs d'une variable numérique. Cependant, le choix de l'intervalle de découpage des valeurs n'est pas toujours aisé. En effet, avec un intervalle trop faible, on fait apparaître trop de variations souvent insignifiantes, tandis qu'avec un intervalle trop élevé les variations de la répartition s'effaçent au profit d'une distribution peu "discriminante" et ne laissent pas apparaître une éventuelle distribution bimodale. Le choix de la taille de l'intervalle s'effectue avec l'option breaks de la fonction hist(). Si cette option n'est pas renseignée, R calcule un intervalle par défaut (option "Sturges") selon laquelle h=range(x)/log2(n)+1. Il existe également d'autres estimations de la taille optimale de l'intervalle de classes à partir d'algorithmes d'optimisation entre biais et variance pour des distributions de référence considérée comme normales [Venables2002] (p. 112). Comme indiqué dans l'aide en ligne (?hist), il faut spécifier la méthode "Scott" ou "FD" (pour "Friedman-Diaconis"). Pour ces méthodes, l'intervalle choisi est calculé comme

h = 3.5*s*n^(-1/3) (Scott, 1979, voir [Venables2002])

ou

h = 2*R*n^(-1/3) (Freedman & Diaconis, 1981, ibidem.).

On peut tester la première méthode à l'aide de données simulées :

# on vérifie l'évolution de h(n)
x <- rnorm(10000)
n <- length(x)
h <- vector("numeric",n-1)
for (i in 2:n) {
  h[i-1] <- (max(x)-min(x))/(log(i,2)+1)
}
plot(h,type="l",xlab="n")
# comparaison entre les choix pour `breaks` et l'option "Sturges"
xx <- seq(-4,4,by=.01)
idx <- c(50,20,8)
par(mfrow=c(2,2))
hist(x,main="Méthode Sturges",ylab="Densité",ylim=c(0,0.4),proba=T)
lines(xx,dnorm(xx),col="red",lwd=2)
legend(1.1,0.4,legend="N(0;1)",lty=1,lwd=2,col="red")
for (i in idx) {
  hist(x,breaks=seq(-4,4,len=i),main=paste("h=",8/i,"(n=50)"),ylim=c(0,0.4),
    ylab="Densité",proba=T)
  lines(xx,dnorm(xx),col="red",lwd=2)
}

Fig. 1 - hist() et l'option breaks

Des estimateurs locaux de densité sont superposés sur chacun des graphiques précédents. Ces estimations de densité (non-paramétrique) sont accessibles à l'aide de la fonction density(), qui accepte pour principaux paramètres :

• la taille de la fenêtre de lissage ;
• le type de noyau à choisir parmi : "gaussian", "epanechnikov", "rectangular","triangular", "biweight", "cosine", et "optcosine".

On peut voir ce que cela donne en variant le paramètre adj :

y <- rnorm(100,0,(1+2*rbinom(100,1,0.35))) # [Venables2002]
par(mfrow=c(1,3))
for (i in seq(0.5,1.5,by=0.5)) {
  hist(y,main=paste("adj=",i),ylim=c(0,0.3),ylab="Densité",proba=T)
  lines(density(y,adj=i),col="blue",lwd=2)
}

Fig. 2 - density() et le paramètre adj

Lorsque le paramètre de lissage est faible (adj=0.5), la courbe de densité souligne beaucoup de variations d'effectifs qui ne semblent pas particulièrement remarquables, tandis qu'à l'inverse, avec un fort lissage (adj=1.5), seule subsiste la tendance générale unimodale. On notera que l'asymétrie gauche reste visible, comme avec adj=1 (option par défaut sous R).

On peut visualiser une démonstration de ce procédé dynamiquement à l'aide du script suivant : density_tk.R. Ce script nécessite la librairie tkrplot.

Voici un autre exemple d'utilisation de l'histogramme sous R :

On peut enfin s'amuser à reproduire un des exemples données par John Fox (ref www ?):

old.par <- par(mar=c(5,4,4,5))
xx <- c(.2, .5, .8)
yy <- xx*0.8+0.1
plot(c(0,1), c(0,1), type="n", axes=FALSE, xlab=expression(Theta), ylab="")
axis(1, at=xx, labels=c(expression(theta[1]),expression(theta[2]),expression(theta[3])))
axis(2, at=yy, las=1, labels=c(expression(widehat(y)[1]),
      expression(widehat(y)[2]),expression(widehat(y)[3])))
axis(4, at=seq(0,1,by=0.5), las=1)
box()
# draw the regression line
abline(0.1,0.8)
# add normal deviates for each Y|x_i plus corresponding probability
# density for a binary response
x <- seq(-3,3,length=100)
y <- dnorm(x)/2
threshold <- 0.35
for (i in 1:3) {
  y1 <- x/12 + yy[i]
  x1 <- y + xx[i]
  lines(x1, y1)
  whichy <- y1>=threshold
  new.x1 <- c(xx[i],x1[whichy])
  new.y1 <- c(y1[whichy][1],y1[whichy])
  polygon(new.x1, new.y1, col="light blue", border=NA)
}
abline(h=yy, lty=3, col="gray")
abline(h=threshold,lwd=3,col="blue")
mtext(c(expression(Upsilon1),expression(Psi)),side=c(2,4),line=3,las=1)
par(old.par)

Fig. 2 - densityplot() (variante)

Graphiques conditionnels

La librairie lattice offre des possibilités graphiques additionnelles puisqu'elle propose les mêmes classes d'outils graphique (histogramme, boîte à moustaches, diagramme, etc.) mais avec la possibilité d'effectuer des représentations graphiques conditionnellement aux valeurs ou niveaux d'une variable concomittante.

library(ade4)
library(lattice)
data(deug)
x <- deug$tab$Algebra
y <- deug$result
densityplot(~x|y,xlab="algèbre",ylab="densité",panel=function(x,...) {
  panel.mathdensity(dmath=dnorm,args=list(mean=mean(x),sd=sd(x)),col="red")
  panel.histogram(x,breaks=NULL,col="cyan")})

Fig. 2 - densityplot()

postscript("densityplot.eps", width = 10.0, height = 10.0,
  horizontal = FALSE, onefile = FALSE, paper = "special",
  encoding = "ISOLatin1.enc")

Fig. x - sexe_scores.png

Diagrammes de dispersion

On peut utiliser ce type de représentation dans le cas où l'on croise plusieurs variables, par exemple:

plot(iris[,1:4],bg=c("red","green3","blue")[iris[,5]],
  pch=c(21,25,24)[iris[,5]],main="Les iris de Fisher",
  labels=c("Longueur\nsépale","Largeur\nsépale",
  "Longueur\npétale","Largeur\npétale"))

Fig. x - Les iris de Fisher

La fonction pairs permet également de représenter une matrice graphique:

a<-rnorm(100)
b<-2*a+rnorm(100)
c<-5*a+rnorm(100)+runif(100)*2
pairs(cbind(a,b,c))

Fig. x - pairs

On peut personnaliser la diagonale de ce graphique grâce à la librairie

Mais dans le cas où le nombre de points est très important (n > 10000), cela devient difficile d'exporter directement l'image à moins d'accepter de gérer des images postscript ou png de 12 Mo !

La libriarie hexbin fournie avec les extensions de Bioconductor (www.bioconductor.org) permet de représenter de larges volumes de données:

library(hexbin)
# adapted from the example in ?hexbin
x <- rnorm(10000)
y <- rnorm(10000)
plot(hexbin(x, y + x*(x+1)/4), ylab="Y",
     main=expression(frac(X %.% (X+1),4) + Y))

Fig. x - hexbin

On retrouve une représentation des effectifs à peu près comparable dans le package survey. Dans un premier temps, on peut représenter les points de coordonnées (x,y) avec une taille qui est proportionnelle à l'effectif associé:

# from ?svyplot
data(api)
dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)

svyplot(api00~api99, design=dstrat, style="bubble")

Fig. x - svyplot

ou alors faire appel à une représentation reposant sur la librarie hexbin grâce à l'option style="hex":

svyplot(api00~api99, design=dstrat, style="hex", xlab="1999 API",
        ylab="2000 API")
svyplot(api00~api99, design=dstrat, style="grayhex",legend=0)

La librarie ggplot (remplacée par ggplot2 à l'heure actuelle) permet également ce type de représentation :

Fig. x - ggplot

En dernier lieu, on peut imager une représentation de densité 2D avec contour ou contourplot (lattice).

Boîtes à moustaches

ns    <- c(15,28,10,20,35)
n     <- length(ns)
group <- factor(rep(1:n,ns),labels=paste("g",1:n,sep=""))
data  <- rnorm(length(group),mean=100+(as.numeric(group)-2)^2)
boxplot(data~group,border=1:n,xlab="Groupe",ylab="Réponse",varwidth=T)
for (i in 1:n) {
  rug(data[as.numeric(group)==i],side=2,col=i)
}

Fig. x - boxplot() + rug()

Il peut être intéressant de superposer sur cette représentation synthétique de la distribution de la variable étudiée la distribution univariée des observations:

x <- rnorm(20,mean=20,sd=2.5)
y <- rnorm(20,mean=22,sd=2.3)

boxplot(x,y)
points(c(rep(1,20),rep(2,20)),c(x,y),col='gray50')
points(c(1,2),c(mean(x),mean(y)),pch='x',cex=2,col=c('blue','red'))

Fig. x - boxplot() + points()

Les paramètres graphiques sont extrêmement faciles à moduler et il est très facile d'obtenir des figures un peu plus "sexy" :

Fig. x - scores_version.png

On peut également utiliser la fonction stripchart qui permet de ne pas superposer les points les uns sur les autres (utile quand il y a beaucou d'observations et que certaines d'entre elles prennent la même valeur), comme ceci:

y <- rnorm(1000)
x <- gl(4,25,1000,labels=paste("g",1:4,sep=""))
stripchart(y~x,method="jitter",ylab="Y",xlab="X",vertical=T)
(moy <- as.numeric(tapply(y,x,mean)))
points(1:4,moy,pch="X",col='blue',cex=2)
lines(1:4,moy,col='blue',lwd=2)

Fig. x - stripchart avec l'option method="jitter"

L'option method="stack" peut également se révéler intéressante:

y <- rpois(500,10)
x <- gl(4,25,500,labels=paste("g",1:4,sep=""))
stripchart(y~x,method="stack",ylab="Y",xlab="X",vertical=T)

Fig. x - stripchart avec l'option method="stack"

Distributions univariées

Les graphiques en barres sont sans doute les plus utilisés dans le cadre de la représentation de la distribution des effectifs en fonction des modalités d'une variable qualitative.

x1 <- c(23.2,34.5,76.3,65.8,12.6)
x2 <- c(15.6,12.4,21.8,20,5.2)
A <- gl(5,1,5,labels=c("a1","a2","a3","a4","a5"))
data <- cbind(x1,x2)
rownames(data) <- levels(A)
barplot(x1,names.arg=levels(A))
barplot(t(data),beside=T,ylim=c(0,100),legend.text=colnames(data),
  col=c("grey50","grey80"),ylab="Fréquence")

Fig. x - barplot()

L'option beside permet de ne pas superposer les différentes catégories, ce qui est parfois plus lisible.

Fig. x - barplot() sans beside=T

On peut également utiliser une boîte à camembert pour représenter les effectifs ventilés sur les modalités de la variable d'intérêt :

names(x1) <- levels(A)
pie(x1/100)

Fig. x - pie()

Remarquons cependant que ces représentations dans lesquelles les proportions relatives sont évaluées par des secteurs angulaires deviennent vite difficiles à analyser lorsqu'il y a beaucoup de modalités (cf. [1]) ou lorsque l'on souhaite croiser différentes variables, et il est préférable d'utiliser des diagrammes en barres.

On peut utiliser cette représentation en barre dans le cas d'une variable continue ventilée sur différentes modalités:

y <- c(rnorm(50,mean=20),rnorm(50,mean=15,sd=3))
x1 <- gl(2,50,100,labels=paste("c",1:2,sep=""))
x2 <- factor(rep(1:4,25),labels=letters[1:4])
y.means <- tapply(y,list(x1,x2),mean)
y.sd <- tapply(y,list(x1,x2),sd)

my.barplot <- function(data, err, ...) {
  tmp <- barplot(data, ...)
  if (length(err) != length(data))
      stop("la longueur de 'data' et 'err' ne correspond pas")
  for (i in 1:length(err))
      arrows(tmp[i],data[i]-err[i],tmp[i],data[i]+err[i],
             code=3,angle=90,length=0.08)
}

my.barplot(y.means,y.sd,beside=T,ylab="Y",xlab="X2",
           ylim=c(0,26),names.arg=levels(x2),legend.text=levels(x1),
           cex.axis=.7,cex.names=.7,cex.lab=.7)

Fig. x - une fonction barplot "améliorée"

On peut encore pousser la personnalisation un peu plus loin:

a <- matrix(NA,nrow=8,ncol=2)
a[,1] <- rev(c(2,3,3,6,5,12,15,20))
a[,2] <- rev(c(1,4,2,3,9,10,30,30))
rownames(a) <- LETTERS[1:8]
colnames(a) <- c("1","2")
cols <- c(rgb(105/255,166/255,233/255),
          rgb(180/255,210/255,244/255),
          rgb(0/255,78/255,162/255))

a.prop <- round(a/sum(a)*100,2)
a.margin <- apply(a,1,sum)
x.offset <- 25

par(mar=c(5,6,4,2))
bp <- barplot(t(a),names.arg=rownames(a),horiz=T,col=cols[1:2],
              las=1,xlab="Effectif",main="",xlim=c(0,max(a.margin)+x.offset),
              border=cols[3])
legend("topright",c("1","2"),pch=rep(19,2),col=cols[1:2],bty="n")
labs <- paste(paste(as.character(a.prop[,1]),"%",sep=""),
              paste(as.character(a.prop[,2]),"%",sep=""),sep=" / ")
text(a.margin,bp,labs,pos=4,cex=.8)

Fig. x - une fonction barplot "améliorée" (2)

ou écrire une fonction qui permettent de représenter conjointement une distribution exprimée en termes d'effectifs et de fréquences:

my.barplot <- function(data, labels, inside=F, signif.digits=1, ...) {
  # Draw a barplot with corresponding frequency (%)
  # * Input *
  #    data   : a table
  #    labels : a vector of labels for y-axis
  #    inside : if TRUE, write the corresponding frequency in the middle
  #             of the bar, else at its end. If the bar is not as large as
  #             the text, then it is written at the end of the bar.
  # * Output *
  #    graphical side-effects
  # * Test *
  #   data <- rbinom(70000,5,0.2); my.barplot(table(data),1:6)
  # 29/11/06, chl
  hoffset <- .025*max(data)
  labs <- round(data/sum(data)*100,signif.digits)
  tmp <- barplot(data, horiz=T, xlim=c(0,max(data)+15*hoffset),
                 names.arg=labels, las=1, ...)
  if (inside) {
    xx <- ifelse(strwidth(labs)<data,data/2-
                 strwidth(paste(labs,"%",sep=" "))/2,data+hoffset)
    text(xx,tmp,paste(labs,"%",sep=" "),pos=4,cex=.8)
  }
  else
    text(data+hoffset,tmp,paste(labs,"%",sep=" "),pos=4,cex=.8)
}

Autre type de personnalisation :

Fig. x - une fonction barplot "améliorée" (3)

Représentations conjointes

Il existe plusieurs fonctions utiles dans le cadre de la représentation des effectifs ventilés sur les modalités de deux variables qualitatives. Ces méthodes de représentation permettent d'éviter l'éventuel superposition des points telle qu'on peut l'observer dans un simple diagramme de dipersion.

Par exemple, la fonction sunflowerplot() représente les effectifs sous forme de diagramme étoilé pour chacun des croisements des modalités des variables :

x <- rpois(500,lambda=2)
y <- rpois(500,lambda=2)
layout(t(matrix (1:2)))
plot(x,y,pch=19)
sunflowerplot(x,y,pch=19)

Fig. x - plot() et sunflowerplot()

Dans le cas de deux variables dichotomiques (e.g. en épidémiologie, exposition vs. maladie), on peut utiliser une représentation en quarts de cercle de l'association entre les deux variables :

x <- c(24.3,75.7,16.8,83.2)
a <- matrix(x,nr=2,byrow=T)
var1 <- c("E-","E+")
var2 <- c("M-","M+")
dimnames(a) <- list(var2,var1)
fourfoldplot(a)

Fig. x - fourfoldplot()

Représentations conditionnelles

Voici un exemple d'utilisation d'un diagramme de fréquence conditionnel (deux variables qualitatives), d'après l'exemple trouvé dans l'aide en ligne (?dotchart) :

dotchart(VADeaths, main = "Death Rates in Virginia - 1940")

Fig. x - dotchart()

La librairie lattice offre également la possibilité de produire des graphiques conditionnels, que le conditionnement porte sur des variables catégorielles ou numériques

Fig. x - score_comp_age.png

Fig. x - score_langue.png