Everyone is welcome here --- except those who have borrowed books from me for and have not returned them yet!

Introduction to R

Posted on May 23, 2015 in stats

Why R (rather than, say, Matlab or Python)?

Pros:

  • advanced statistics (i.e. hierarchical mixed-effects modeling)
  • powerful data manipulation (replaced awk/sed/grep for me)
  • help reproducible science (scripting, free, multiplatform software)
  • easy to generate complex data for simulations
  • very nice reporting with Seave or rmardown (like ipython notebooks but better IMHO)
  • lots of documentation available
  • high-level symbolic objects (compared to Matlab which is only good with numbers)

Cons

  • it requires learning yet another programming language
  • the language is a bit exotic and rich (R is a clone of S invented in 1998), but one can use only a subset of the language.
  • The online documentation is cryptic (like unix man pages)

Disclaimer: I have not used nor Python statsmodels (see Edouard Duchesnay), nor the new Matlab statistical toolbox (the old one was a scandal),

Installation du système de base

Les instructions de téléchargement et d’installation sont disponibles sur le site principal du logiciel R (http://www.r-project.org).

Je recommande vivement d’installer l’environement de travail rstudio (http://www.rstudio.com).

Les inconditionnels d’Emacs apprécieront les packages ESS et polymode.

On m’a dit du bien de https://jasp-stats.org/.

Modules additionnels

De nombreux packages ajoutant des fonctions à R sont disponibles sur Contributed extension packages du site CRAN. Par exemple:

  • ez: Easy Analysis and Visualization of Factorial Experiments
  • nlme, lme4: Linear Mixed-Effects Models
  • vcd : Visualizing Categorical Data
  • psy : Various procedures used in psychometry
  • multcomp: Corrections for multiple comparisons
  • ggplot2: Elegant Data Visualisations Using the Grammar of Graphics
  • plyr, reshape: Tools for Splitting, Applying and Combining Data

Ils peuvent être installés avec install.packages() dans R, par des menus, ou bien, sous Linux, directement en ligne de commande R CMD INSTALL package.tar.gz.

Documentation

  • Manuels officiels de R “An Introduction to R” et “R Data Import et Export” (ok for people who already know how to program in other languages)

  • Documentation/Other/Contributed" on the CRAN (see https://cran.r-project.org/other-docs.html). For example:

  • R pour les débutants par Emmanuel Paradis.
  • Introduction au système R par Yves Brostaux.

  • le site Statistiques avec R, réalisé par Vincent Zoonekynd, à l’adresse suivante :
    http://zoonek2.free.fr/UNIX/48_R/all.html.

  • Books

  • Introductory statistics with R par Peter Dalgaard
  • Psychologie statistique avec R par Yvonnick Noël
  • Modeling Psychophysical Data in R par Knoblauch & Maloney
  • An R and S-plus companion to applied regression par John Fox
  • Regression Modeling Strategies par Franck Harrell
  • Analyzing Linguistic Data par Harald Baayen
  • Data Analysis Using Regression and Multilevel / Hierarchical Models par Gelman & Hill

Premiers pas

Ouvrir rstudio et sélectionner ‘New/Rmarkdown’ dans le menu ‘File’.

Un document rmarkdown (Rmd) mélange du texte (au format markdown) et du code R, placé entre à l’intérieur de chunks, délimités par des backtics.

Le button knit génère un rapport au format html, pdf ou Word.

Dans Rstudio, pour insérer du code R dans un document rmarkdown, utiliser le menu ‘Code/Insert chunk’ (ou ‘Ctrl-Alt-I’)

```r
2 + 3
```

```
## [1] 5
```

Poursuivez avec:

a = 1:10
a
##  [1]  1  2  3  4  5  6  7  8  9 10
b = rnorm(10)
plot(a, b, pch=16, col=2, type='b')

a=c(3, 4, 6, 7, 8, 9) # vecteur numérique à six éléments.
a
## [1] 3 4 6 7 8 9
length(a)
## [1] 6
b = c('alpha', 'beta') # vecteur contenant deux chaînes de caractères.
b
## [1] "alpha" "beta"
length(b)
## [1] 2

La liste des variables dans le workspace courant peut être affichée par ls(), et une variable peut être détruite par la commande rm(nom).

Entrez les commandes suivantes, pas à pas, et observez le résultat:

a = rnorm(20, mean=55, sd=10)
mean(a)
## [1] 59.62146
sd(a) 
## [1] 10.48517
max(a)
## [1] 75.89468
summary(a)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   41.05   51.86   62.15   59.62   67.59   75.89
hist(a)

boxplot(a)

stripchart(a, pch=16, cex=2, col=2, method='jitter', vertical=T)

x1 = rnorm(30, mean=100, sd=10)
x2 = rnorm(30, mean=110, sd=10)
boxplot(x1, x2)

t.test(x1, x2)
## 
##  Welch Two Sample t-test
## 
## data:  x1 and x2
## t = -4.027, df = 57.837, p-value = 0.0001666
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -13.989379  -4.699221
## sample estimates:
## mean of x mean of y 
##  100.2899  109.6342
y = x1 + 4 * x2 + rnorm(30, sd=30)
plot(y ~ x2)
summary(a <- lm(y ~  x2))
## 
## Call:
## lm(formula = y ~ x2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -69.982 -18.354  -4.733  20.051  63.009 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  50.6459    73.2149   0.692    0.495    
## x2            4.4198     0.6658   6.639 3.34e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.35 on 28 degrees of freedom
## Multiple R-squared:  0.6115, Adjusted R-squared:  0.5976 
## F-statistic: 44.07 on 1 and 28 DF,  p-value: 3.342e-07
abline(a)

Aide en ligne

A tout moment, une aide en ligne est disponible à l’aide de la commande help.search(‘mot clé’). La description détaillée d’une commande s’obtient en tapant ?nom de la commande’.

Essayez :

?t.test
apropos('t.test')
help.search("test")

Quitter R

Quitter R et répondez “Oui” à la question “Save workspace image?”

Redémarrez R, et remarquez la ligne:

[Previously saved workspace restored]

Tapez “ls()” et constatez que vos variables sont toujours là.

Le “workspace” (“espace de travail”), c’est à dire l’ensemble des variables, a été sauvegardé sur le disque. Cela permet de reprendre une analyse de données au point où on l’a laissée quand on a quitté R.

Si vous voulez “nettoyer” le workspace, c’est à dire supprimer toutes les variables qu’il contient, tapez la commande “rm(list=ls())”.

Dans Rstudio, utilisez les projets! (un répertoire par projet avec données et code; permet de passer d’une projet à l’autre sans interférence)

Manipulations de base

Vecteurs

L’objet de base en R est le vecteur. Un vecteur peut contenir des valeurs numériques, des valeurs de vérité (True or False), des chaînes de caractères… Les fonctions les plus utilisées pour créer des vecteurs sont c, rep et seq :

c(1, 2, 3, 4, 5, 6) 
## [1] 1 2 3 4 5 6
c(T, T, F, F)
## [1]  TRUE  TRUE FALSE FALSE
c('a', 'b')
## [1] "a" "b"
rep(55, 10)
##  [1] 55 55 55 55 55 55 55 55 55 55
rep(c(1,2), 10)
##  [1] 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2
rep(c('a', 'b'), c(2, 7))
## [1] "a" "a" "b" "b" "b" "b" "b" "b" "b"
seq(1, 10, by=.1)
##  [1]  1.0  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9  2.0  2.1  2.2  2.3
## [15]  2.4  2.5  2.6  2.7  2.8  2.9  3.0  3.1  3.2  3.3  3.4  3.5  3.6  3.7
## [29]  3.8  3.9  4.0  4.1  4.2  4.3  4.4  4.5  4.6  4.7  4.8  4.9  5.0  5.1
## [43]  5.2  5.3  5.4  5.5  5.6  5.7  5.8  5.9  6.0  6.1  6.2  6.3  6.4  6.5
## [57]  6.6  6.7  6.8  6.9  7.0  7.1  7.2  7.3  7.4  7.5  7.6  7.7  7.8  7.9
## [71]  8.0  8.1  8.2  8.3  8.4  8.5  8.6  8.7  8.8  8.9  9.0  9.1  9.2  9.3
## [85]  9.4  9.5  9.6  9.7  9.8  9.9 10.0

Factors

Les facteurs sont un type de vecteur utilisés pour classifier les valeurs d’autres vecteurs (les facteurs sont des “variables indicatrices” de sous-groupe).

On peut créer un facteur à partir d’un vecteur grâce à la fonction factor, ou directement avec la fonction “gl”.

(a=factor(rep(c('alpha', 'beta'), c(10, 10))))
##  [1] alpha alpha alpha alpha alpha alpha alpha alpha alpha alpha beta 
## [12] beta  beta  beta  beta  beta  beta  beta  beta  beta 
## Levels: alpha beta
(b=gl(3, 4, 48, labels=c('a', 'b', 'c')))
##  [1] a a a a b b b b c c c c a a a a b b b b c c c c a a a a b b b b c c c
## [36] c a a a a b b b b c c c c
## Levels: a b c
(x=rnorm(48))
##  [1]  0.08663439  0.56640818 -1.28564389 -0.46409367  0.61855196
##  [6] -1.73643980 -1.01154095  0.81757449 -0.19542090 -1.70557052
## [11] -0.71099622  0.28185372  0.73250364  0.45536062 -0.65901928
## [16]  1.53296209 -0.83867281 -0.23856889  0.79166088 -0.31387160
## [21]  0.59093975 -0.09839616 -0.85343102 -0.61579513 -0.39374777
## [26]  3.03846985 -0.57356106  1.19269667 -0.11507413  1.66939049
## [31] -0.11033377 -0.78941935 -0.26974258 -0.41992362  0.09838458
## [36]  2.14619194  0.32666938  1.31231207  0.34092964 -2.76864549
## [41] -0.80538602 -0.29854732  1.06787524 -0.50644784 -0.56554538
## [46] -0.25865954 -0.19724045 -1.29884523
tapply(x, b, mean)
##          a          b          c 
##  0.2150147 -0.1124531 -0.2545123
boxplot(x ~ b)

stripchart(x ~ b, method='jitter')

stripchart(x ~ b, method='jitter', vertical=T)

Accéder aux éléments d’un vecteur

(a = rnorm(50))
##  [1]  1.25921678 -0.23884975 -0.13484160 -1.51478143 -0.56919971
##  [6]  0.91664171  0.64844004  0.78537044 -0.18665961  0.78144962
## [11]  0.17210125 -1.20155180 -0.53363569  0.07508076  0.76674842
## [16] -0.51109724 -2.55442168 -1.47833924 -0.55419893  0.52625647
## [21]  0.80573352  0.07994848  0.48957025 -1.24791692  1.22014045
## [26] -0.05253264 -0.59210435  0.19555335 -1.46052462  0.07783316
## [31] -0.63070639  0.31288194 -0.27525230  0.98934769 -0.83328276
## [36] -0.92679672  1.12147618 -0.49450154  1.02259277  0.01522365
## [41] -0.30102486  1.41664846 -0.17527009  0.85976801 -0.13185991
## [46]  1.28885322 -0.63915249 -1.72396133 -0.35864218 -1.49272223
a[1]
## [1] 1.259217
a[2]
## [1] -0.2388498
a[c(1, 3, 5)]
## [1]  1.2592168 -0.1348416 -0.5691997
a>0
##  [1]  TRUE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE
## [12] FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE
## [23]  TRUE FALSE  TRUE FALSE FALSE  TRUE FALSE  TRUE FALSE  TRUE FALSE
## [34]  TRUE FALSE FALSE  TRUE FALSE  TRUE  TRUE FALSE  TRUE FALSE  TRUE
## [45] FALSE  TRUE FALSE FALSE FALSE FALSE
a[a>0]
##  [1] 1.25921678 0.91664171 0.64844004 0.78537044 0.78144962 0.17210125
##  [7] 0.07508076 0.76674842 0.52625647 0.80573352 0.07994848 0.48957025
## [13] 1.22014045 0.19555335 0.07783316 0.31288194 0.98934769 1.12147618
## [19] 1.02259277 0.01522365 1.41664846 0.85976801 1.28885322
a[abs(a-mean(a))<2*sd(a)]
##  [1]  1.25921678 -0.23884975 -0.13484160 -1.51478143 -0.56919971
##  [6]  0.91664171  0.64844004  0.78537044 -0.18665961  0.78144962
## [11]  0.17210125 -1.20155180 -0.53363569  0.07508076  0.76674842
## [16] -0.51109724 -1.47833924 -0.55419893  0.52625647  0.80573352
## [21]  0.07994848  0.48957025 -1.24791692  1.22014045 -0.05253264
## [26] -0.59210435  0.19555335 -1.46052462  0.07783316 -0.63070639
## [31]  0.31288194 -0.27525230  0.98934769 -0.83328276 -0.92679672
## [36]  1.12147618 -0.49450154  1.02259277  0.01522365 -0.30102486
## [41]  1.41664846 -0.17527009  0.85976801 -0.13185991  1.28885322
## [46] -0.63915249 -1.72396133 -0.35864218 -1.49272223
(b=gl(2, 25, labels=c('g1', 'g2')))
##  [1] g1 g1 g1 g1 g1 g1 g1 g1 g1 g1 g1 g1 g1 g1 g1 g1 g1 g1 g1 g1 g1 g1 g1
## [24] g1 g1 g2 g2 g2 g2 g2 g2 g2 g2 g2 g2 g2 g2 g2 g2 g2 g2 g2 g2 g2 g2 g2
## [47] g2 g2 g2 g2
## Levels: g1 g2
a[b=='g1']
##  [1]  1.25921678 -0.23884975 -0.13484160 -1.51478143 -0.56919971
##  [6]  0.91664171  0.64844004  0.78537044 -0.18665961  0.78144962
## [11]  0.17210125 -1.20155180 -0.53363569  0.07508076  0.76674842
## [16] -0.51109724 -2.55442168 -1.47833924 -0.55419893  0.52625647
## [21]  0.80573352  0.07994848  0.48957025 -1.24791692  1.22014045

Une particularité de R est que les éléments d’un vecteur peuvent avoir des noms:

v=c(1, 2, 3, 4)
names(v) = c('alpha', 'beta', 'gamma', 'delta')
v['beta']
## beta 
##    2

Cela permet de créer des dictionnaires. Par exemple, un vecteur ‘freq’ donnant la fréquence d’usage des mots peut avoir les mots comme ‘names’ ; il suffit alors de taper “freq['aller']” pour obtenir la fréquence du mot ‘aller’.

freq=c(45, 3)
names(freq) = c('aller', 'vaquer')
freq
##  aller vaquer 
##     45      3
freq['aller']
## aller 
##    45

Arrays, listes et data.frames

D’autres objets de R sont: - les listes (comme un cell array en, Matlab, mais avec des noms pour les éléments) - les arrays (comme numpy.array en Python ou array en Matlab) - les data.frames (liste de vecteurs de même longueur represetant des tables de données)

(a=array(1:20, dim=c(4, 5)))
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    5    9   13   17
## [2,]    2    6   10   14   18
## [3,]    3    7   11   15   19
## [4,]    4    8   12   16   20
a[2, 4]
## [1] 14
(b=list(alpha=1:3,
        beta=c('a','b','c','d')))
## $alpha
## [1] 1 2 3
## 
## $beta
## [1] "a" "b" "c" "d"
names(b)
## [1] "alpha" "beta"
b$alpha
## [1] 1 2 3
b$beta
## [1] "a" "b" "c" "d"
(c=data.frame(a=gl(2, 5, 10),
              b=1:10,
              x=rnorm(10)))
##    a  b           x
## 1  1  1  0.04143612
## 2  1  2  0.45721924
## 3  1  3  1.23314302
## 4  1  4 -1.55660935
## 5  1  5 -0.42638813
## 6  2  6  0.25028516
## 7  2  7 -0.25124132
## 8  2  8  0.71585284
## 9  2  9 -0.76836867
## 10 2 10 -0.63521548
c$a
##  [1] 1 1 1 1 1 2 2 2 2 2
## Levels: 1 2
c$b
##  [1]  1  2  3  4  5  6  7  8  9 10
c$x
##  [1]  0.04143612  0.45721924  1.23314302 -1.55660935 -0.42638813
##  [6]  0.25028516 -0.25124132  0.71585284 -0.76836867 -0.63521548
c[1:2,]
##   a b          x
## 1 1 1 0.04143612
## 2 1 2 0.45721924

Deux fonctions essentielles: tapply et aggregate

dat = data.frame(a = gl(2, 10, 40),
                 b = gl(4, 5, 40),
                 y = rnorm(40))

with(dat, tapply(y, list(a=a, b=b), mean))
##    b
## a            1           2           3          4
##   1 -0.3927085 0.006901131          NA         NA
##   2         NA          NA -0.07677625 -0.0441236
with(dat, aggregate(y, list(a=a, b=b), mean))
##   a b            x
## 1 1 1 -0.392708498
## 2 1 2  0.006901131
## 3 2 3 -0.076776246
## 4 2 4 -0.044123597

Variables

Les objets peuvent être enregistrés dans des variables avec l’opérateur = (ou <-). Pour voir le contenu de l’objet représenté par une variable, il suffit de taper le nom de celle-ci.

a<-c(1,2,3)
a
## [1] 1 2 3
ls()
##  [1] "a"    "b"    "c"    "dat"  "freq" "v"    "x"    "x1"   "x2"   "y"
rm(a)
ls()
## [1] "b"    "c"    "dat"  "freq" "v"    "x"    "x1"   "x2"   "y"

Les vecteurs contenus dans une liste ou dans un data.frame sont accessibles avec le symbole $. Un data.frame peut être “attaché” pour que ses vecteurs soient directement accessibles.

mydata <- data.frame(a=gl(2, 5, 10),
                     b=1:10,
                     x=rnorm(10))
names(mydata)
## [1] "a" "b" "x"
mydata$a
##  [1] 1 1 1 1 1 2 2 2 2 2
## Levels: 1 2
mydata$b
##  [1]  1  2  3  4  5  6  7  8  9 10
mydata$x
##  [1] -0.13643701  0.06632392  1.53001350 -0.91350929 -0.29701974
##  [6] -1.46132006 -0.66631967 -1.43775172 -0.11631767 -1.77773141
attach(mydata)
## The following objects are masked _by_ .GlobalEnv:
## 
##     b, x
a
##  [1] 1 1 1 1 1 2 2 2 2 2
## Levels: 1 2
b
## $alpha
## [1] 1 2 3
## 
## $beta
## [1] "a" "b" "c" "d"
x
##  [1]  0.08663439  0.56640818 -1.28564389 -0.46409367  0.61855196
##  [6] -1.73643980 -1.01154095  0.81757449 -0.19542090 -1.70557052
## [11] -0.71099622  0.28185372  0.73250364  0.45536062 -0.65901928
## [16]  1.53296209 -0.83867281 -0.23856889  0.79166088 -0.31387160
## [21]  0.59093975 -0.09839616 -0.85343102 -0.61579513 -0.39374777
## [26]  3.03846985 -0.57356106  1.19269667 -0.11507413  1.66939049
## [31] -0.11033377 -0.78941935 -0.26974258 -0.41992362  0.09838458
## [36]  2.14619194  0.32666938  1.31231207  0.34092964 -2.76864549
## [41] -0.80538602 -0.29854732  1.06787524 -0.50644784 -0.56554538
## [46] -0.25865954 -0.19724045 -1.29884523
detach(mydata)

Lire des données

Quand les données sont très peu nombreuses, on peut les entrer directement dans un vecteur (comme on l’a fait jusqu’ici) avec la fonction ‘c’.

Les fonctions scan, read.tableet read.csv permettent de lire des données enregistrées dans des fichiers textes.

scan lit une suite de données dans un vecteur.

Avec un éditeur de texte, créez un fichier datafile1.txt contenant:

3.4 5.6 2.1 6.7 8.9

Puis, dans R, entrez:

scores<-scan('datafile1.txt')

On peut également entrer des données directement en ligne de commande :

scores<-scan('')

La fonction read.table lit des données présentées sous forme tabulaire (par ex. les fichiers .csv enregistrés par Excel) et renvoit un data.frame.

Créez un fichier datafile2.txt contenant:

sujet groupe score
s1    exp    3
s2    exp    4
s3    exp    6
s4    cont   7
s5    cont   8

Puis importez le dans R:

a<-read.table('datafile2.txt',header=T)
a

scan et read.table ne lisent que des fichiers textes, Le package ‘foreign’ permet de lire directement certains fichiers de données binaires provenant de SPSS, SAS, …

library(help='foreign')

Mentionnons également l’existence de packages permettant d’accéder à des informations stockées dans des bases de données (MySQL, Oracle…).

Statistiques élémentaires

Cette section a pour but d’illustrer quelques concepts fondamentaux de la statistique inférentielle, et de présenter les principales fonctions de R pour le traitement statistique des données recueillies lors d’un protocole expérimental.

Manipulation des distributions de probabilités

Distributions univariées

Différentes fonctions permettent de générer des nombres aléatoires suivant une certaine distribution de probabilité :

runif(10) # distribution uniforme
##  [1] 0.28955420 0.99407582 0.97490284 0.26175089 0.19499848 0.09255223
##  [7] 0.71100801 0.47537991 0.58177741 0.28756436
rnorm(10) # distribution normale
##  [1] -0.22407969  0.53923472 -0.81184444 -0.64999117 -0.00258321
##  [6] -0.42897630  0.72392279  0.91105106 -1.22639625 -0.70894392
rnorm(10, mean=100)
##  [1]  99.45777 100.78720 100.10059 100.10119 101.39015  98.38287 100.40605
##  [8]  98.76283 100.66653 100.20743
rbinom(10, size=1, prob=.5) # distribution binomiale
##  [1] 1 1 0 0 0 0 1 0 0 0

La fonction rnorm génère des nombres aléatoires distribués selon une loi normale. En augmentant le nombre d’échantillons générés (de 10 à 10000), on constate que la distribution des valeurs obtenues se rapproche de plus en plus d’une distribution normale continue :

s1 = rnorm(10, mean=2)
summary(s1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.5073  1.3550  1.5660  1.7850  2.3190  3.0200
s2 = rnorm(100, mean=2)
summary(s2)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.04981 1.51700 2.18100 2.12500 2.89700 4.70600
s3 = rnorm(10000, mean=2)
summary(s3)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -1.423   1.312   1.989   1.989   2.658   5.572
par(mfrow=c(3,3))       # organisation des graphiques selon une matrice 3 x 3

hist(s1)                # histogrammes
hist(s2)
hist(s3)

                                        # graphes en évitant le chevauchement des points de même coordonnées
stripchart(s1,method='jitter',vert=T,pch=16)    
stripchart(s2,method='jitter',vert=T,pch=16)
stripchart(s3,method='jitter',vert=T,pch='.')


plot(density(s1))       # fonction de densité
x=seq(-5, 5, by=.01)      # vecteur de coordonnées normées pour les abscisses
lines(x, dnorm(x,mean=2), col=2)
plot(density(s2))
lines(x, dnorm(x,mean=2), col=2)
plot(density(s3))
lines(x, dnorm(x,mean=2), col=2)

En première approximation, la distribution théorique de la taille des individus de sexe masculin, français, et dans la tranche d’âge 20-35 ans,suit une loi normale de moyenne 170 et d’écart-type 10.

On peut donc non seulement situer un individu, ou un groupe d’individus, dans cette distribution, mais également évaluer la probabilité qu’un individu choisi au hasard parmi la population entière mesure moins de 185 cm, ou plus de 198 cm, ou ait une taille comprise entre 174 et 186 cm.

Lorsque l’on ne dispose pas des tables de lois normales N(m;s2) (il y en a une infinité puisqu’il y a 2 paramètres libres), on utilise la loi normale centrée-réduite N(0;12) (encore appelée loi Z), dont la table est disponible la plupart des manuels ou bien sur le web. Cependant R fournit directement les tables des lois normales, par l’intermédiaire de la commande pnorm, qui prend en arguments la valeur repère, la moyenne et l’écart-type théoriques.

taille=seq(130, 210, by=1)
plot(taille, dnorm(taille, mean=170, sd=10), type='b', col="red")
p = pnorm(185, mean=170, sd=10)
abline(v=185, col=4)
text(185, .012, paste("P(X<185)=", signif(p, 3)), col=4, pos=2, cex=.6)
p = pnorm(198, mean=170, sd=10)
abline(v=198, col=4)
text(198, .002, paste("P(X>198)=",round(1 - p, 3)), col=4, pos=4, cex=.6)

La probabilité qu’un individu choisi au hasard parmi la population entière mesure moins de 185 cm ( P(X < 185) ) est de 0.933 (obtenu par pnorm(185, mean=170, sd=10)). La probabilité qu’un individu mesure plus de 198 cm est de 0.003 (1-P(X < 198 ), et la probabilité que sa taille soit comprise entre 174 et 186 est 0.290 ( P(X > 186)-P(X < 174) ).

On constate que la probabilité qu’un individu choisi aléatoirement dans une population de moyenne 170 ± 10 mesure plus de 198 cm est très faible. C’est sur la base de ce calcul de probabilités que repose le test de typicalité, ou “test Z” : un groupe d’individus (i.e. un échantillon) sera déclaré atypique ou non représentatif de la population parente dont il est issu, lorsqu’il a une position au moins aussi extrême qu’une certaine position de référence, correspondant en général à la probabilité 0.05.

R permet également de générer d’autres distributions de probabilités, notamment la loi binomiale, les lois statistiques telles que le t de Student, le F de Fisher-Snedecor, le chi-deux ( X^2 ), etc. On peut ainsi voir dans l’exemple qui suit que la distribution du t de Student tend vers la loi normale lorsque la taille de l’échantillon est suffisamment grande (dans cet exemple, on a manipulé le degré de liberté df, donné en argument de la fonction dt).

?pnorm
?pt
?pbinom
help.search('distribution')
pnorm(2)
pt(3,df=10)             # fonction de répartition de la loi du t de Student
qnorm(.99)              # donne la valeur associée au 99ème centile d'une distribution normale
t <- -50:50/10
plot(dnorm(t),type='l',col='red')
par(new=T)              # le prochain graphe sera superposé au précédent
plot(dt(t,df=5),type='l')

Imaginons que vous disposiez d’une pièce dont vous vous demandez si elle est baisée. Vous prévoyez de la lancer 10 fois à pile ou face. A partir de quelle proportion relative d’essais face/pile (ou l’inverse) considérerez- vous que la pièces est truquée ?

Si la pièce n’est pas truquée, le nombre de “pile” suit une loi binomiale.

plot(dbinom(0:10, rep(10,11), prob=1/2), type='h')

hist(rbinom(100, 10, .5))

hist(rbinom(1000, 10, .5))

hist(rbinom(10000, 10, .5))

Supposez que vous tiriez à pile ou face 10 fois de suite, et que la pièce retombe 8 fois sur ‘pile’. Quelle la probabilité d’observer cela si la pièce n’est pas biaisée ?

binom.test(8,10)
## 
##  Exact binomial test
## 
## data:  8 and 10
## number of successes = 8, number of trials = 10, p-value = 0.1094
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.4439045 0.9747893
## sample estimates:
## probability of success 
##                    0.8
prop.test(8,10,1/2) # test approché
## 
##  1-sample proportions test with continuity correction
## 
## data:  8 out of 10, null probability 1/2
## X-squared = 2.5, df = 1, p-value = 0.1138
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.4421814 0.9645731
## sample estimates:
##   p 
## 0.8

Distributions conjointes

Si l’on reprend l’exemple précédent des tailles de la population française masculine (20-25 ans), on a une distribution similaire (i.e. suivant une loi normale de moyenne 70 et d’écart-type 7) pour les poids. On peut bien évidemment se poser les mêmes questions que précédemment, mais on peut également s’intéresser à la relation entre ces deux variables quantitatives. En représentant le poids en fonction de la taille, on peut évaluer la liaison linéaire entre ces deux variables à l’aide du coefficient de corrélation de Bravais-Pearson.

Pour illustrer cela, nous allons utiliser les données issues d’une population d’enfants de sexe masculin âgés de 11 à 16 ans.

taille <- c(172, 155, 160, 142, 157, 142, 148, 180, 167, 165)
poids <-  c(50.5, 38.1, 57.3, 39.3, 46.1, 37.1, 45.9, 66.3, 60, 50.5)
plot(poids ~ taille)
r <- lm(poids ~ taille)     # modèle linéaire (x,y)
summary(r)                      # diagnostic de la régression
## 
## Call:
## lm(formula = poids ~ taille)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.5140 -2.4687  0.1249  3.7292  7.4018 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -55.1964    23.5120  -2.348  0.04686 * 
## taille        0.6568     0.1476   4.449  0.00214 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.582 on 8 degrees of freedom
## Multiple R-squared:  0.7122, Adjusted R-squared:  0.6762 
## F-statistic: 19.79 on 1 and 8 DF,  p-value: 0.002143
abline(r)                       # tracé de la droite de régression

predict(r, list(taille=c(175)))
##        1 
## 59.75083

Ensuite, à partir de la connaissance de cette liaison linéaire, on peut se demander quelle serait le poids théorique (non observé) d’un individu dont on ne connaît que la taille : c’est le domaine de la régression linéaire. L’affichage des paramètres de la droite de régression donne la relation poids = 0.657 x taille - 55.196. Ainsi, on peut prédire que le poids d’un enfant mesurant 175 cm sera de 59.8 kg.

Résumés numériques et représentations graphiques

Résumés numériques

Le résumé statistique des principaux indicateurs descriptifs de position et de dispersion peut être obtenu à l’aide des fonctions mean, sd, median ; la fonction summary donne un résumé plus complet - par exemple, lorsqu’il s’agit d’un vecteur, elle indique la moyenne et la médiane, ainsi que l’étendue et les valeurs des premier et troisième quartiles.

a<-rnorm(100)
mean(a)
## [1] -0.08327925
sd(a) # écart-type corrigé
## [1] 0.9713326
summary(a)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -2.986000 -0.658000 -0.005623 -0.083280  0.475000  2.262000
boxplot(a)

mean(a, trim=.1) # moyenne sans les 10 % d'observations en fin de vecteur
## [1] -0.06920461

Représentations graphiques

Les fonctions graphiques standard en 2D - boxplot, plot, hist - ont été vues dans les sections précédentes. La création de graphiques personnalisés sous R est facilitée par son extrême souplesse quant au paramétrage des graphiques (positionnement, symboles et type de tracés, etc.). L’utilisation de l’aide en ligne est vivement recommandée.

Pour les graphiques en trois dimensions (z étant une matrice de dim 3), on pourra utiliser les fonctions image et contour :

x = 1:10
y = 1:10
z = outer(x, y, "*")
persp(x, y, z)

image(z)

contour(z)

Définition de fonctions

Il est possible de définir ses propres fonctions sous Ret d’enrichir ainsi le langage.

Par exemple, R ne possède pas de fonction pour calculer l’erreur-type. On peut en définir une de la manière suivante :

se <- function (x) { sd(x)/sqrt(length(x)) }

L’exemple suivant permet de calculer la moyenne arithmétique après suppression des valeurs atypiques, i.e. supérieures à 2 écart-types de la moyenne :

clmean <- function (x) {
      m<-mean(x)
      d<-sqrt(var(x))
      threshold<-2
      mean(x[(x-m)/d<threshold])
    }

a<-c(rnorm(100),5)
mean(a)
## [1] -0.02242767
clmean(a)
## [1] -0.07265195

On peut lire le code des fonctions existantes:

clmean
## function (x) {
##       m<-mean(x)
##       d<-sqrt(var(x))
##       threshold<-2
##       mean(x[(x-m)/d<threshold])
##     }
ls
## function (name, pos = -1L, envir = as.environment(pos), all.names = FALSE, 
##     pattern, sorted = TRUE) 
## {
##     if (!missing(name)) {
##         pos <- tryCatch(name, error = function(e) e)
##         if (inherits(pos, "error")) {
##             name <- substitute(name)
##             if (!is.character(name)) 
##                 name <- deparse(name)
##             warning(gettextf("%s converted to character string", 
##                 sQuote(name)), domain = NA)
##             pos <- name
##         }
##     }
##     all.names <- .Internal(ls(envir, all.names, sorted))
##     if (!missing(pattern)) {
##         if ((ll <- length(grep("[", pattern, fixed = TRUE))) && 
##             ll != length(grep("]", pattern, fixed = TRUE))) {
##             if (pattern == "[") {
##                 pattern <- "\\["
##                 warning("replaced regular expression pattern '[' by  '\\\\['")
##             }
##             else if (length(grep("[^\\\\]\\[<-", pattern))) {
##                 pattern <- sub("\\[<-", "\\\\\\[<-", pattern)
##                 warning("replaced '[<-' by '\\\\[<-' in regular expression pattern")
##             }
##         }
##         grep(pattern, all.names, value = TRUE)
##     }
##     else all.names
## }
## <bytecode: 0x2ddbd78>
## <environment: namespace:base>
t.test
## function (x, ...) 
## UseMethod("t.test")
## <bytecode: 0x3639860>
## <environment: namespace:stats>
methods(t.test)
## [1] t.test.default* t.test.formula*
## see '?methods' for accessing help and source code
getAnywhere(t.test.default)
## A single object matching 't.test.default' was found
## It was found in the following places
##   registered S3 method for t.test from namespace stats
##   namespace:stats
## with value
## 
## function (x, y = NULL, alternative = c("two.sided", "less", "greater"), 
##     mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95, 
##     ...) 
## {
##     alternative <- match.arg(alternative)
##     if (!missing(mu) && (length(mu) != 1 || is.na(mu))) 
##         stop("'mu' must be a single number")
##     if (!missing(conf.level) && (length(conf.level) != 1 || !is.finite(conf.level) || 
##         conf.level < 0 || conf.level > 1)) 
##         stop("'conf.level' must be a single number between 0 and 1")
##     if (!is.null(y)) {
##         dname <- paste(deparse(substitute(x)), "and", deparse(substitute(y)))
##         if (paired) 
##             xok <- yok <- complete.cases(x, y)
##         else {
##             yok <- !is.na(y)
##             xok <- !is.na(x)
##         }
##         y <- y[yok]
##     }
##     else {
##         dname <- deparse(substitute(x))
##         if (paired) 
##             stop("'y' is missing for paired test")
##         xok <- !is.na(x)
##         yok <- NULL
##     }
##     x <- x[xok]
##     if (paired) {
##         x <- x - y
##         y <- NULL
##     }
##     nx <- length(x)
##     mx <- mean(x)
##     vx <- var(x)
##     if (is.null(y)) {
##         if (nx < 2) 
##             stop("not enough 'x' observations")
##         df <- nx - 1
##         stderr <- sqrt(vx/nx)
##         if (stderr < 10 * .Machine$double.eps * abs(mx)) 
##             stop("data are essentially constant")
##         tstat <- (mx - mu)/stderr
##         method <- if (paired) 
##             "Paired t-test"
##         else "One Sample t-test"
##         estimate <- setNames(mx, if (paired) 
##             "mean of the differences"
##         else "mean of x")
##     }
##     else {
##         ny <- length(y)
##         if (nx < 1 || (!var.equal && nx < 2)) 
##             stop("not enough 'x' observations")
##         if (ny < 1 || (!var.equal && ny < 2)) 
##             stop("not enough 'y' observations")
##         if (var.equal && nx + ny < 3) 
##             stop("not enough observations")
##         my <- mean(y)
##         vy <- var(y)
##         method <- paste(if (!var.equal) 
##             "Welch", "Two Sample t-test")
##         estimate <- c(mx, my)
##         names(estimate) <- c("mean of x", "mean of y")
##         if (var.equal) {
##             df <- nx + ny - 2
##             v <- 0
##             if (nx > 1) 
##                 v <- v + (nx - 1) * vx
##             if (ny > 1) 
##                 v <- v + (ny - 1) * vy
##             v <- v/df
##             stderr <- sqrt(v * (1/nx + 1/ny))
##         }
##         else {
##             stderrx <- sqrt(vx/nx)
##             stderry <- sqrt(vy/ny)
##             stderr <- sqrt(stderrx^2 + stderry^2)
##             df <- stderr^4/(stderrx^4/(nx - 1) + stderry^4/(ny - 
##                 1))
##         }
##         if (stderr < 10 * .Machine$double.eps * max(abs(mx), 
##             abs(my))) 
##             stop("data are essentially constant")
##         tstat <- (mx - my - mu)/stderr
##     }
##     if (alternative == "less") {
##         pval <- pt(tstat, df)
##         cint <- c(-Inf, tstat + qt(conf.level, df))
##     }
##     else if (alternative == "greater") {
##         pval <- pt(tstat, df, lower.tail = FALSE)
##         cint <- c(tstat - qt(conf.level, df), Inf)
##     }
##     else {
##         pval <- 2 * pt(-abs(tstat), df)
##         alpha <- 1 - conf.level
##         cint <- qt(1 - alpha/2, df)
##         cint <- tstat + c(-cint, cint)
##     }
##     cint <- mu + cint * stderr
##     names(tstat) <- "t"
##     names(df) <- "df"
##     names(mu) <- if (paired || !is.null(y)) 
##         "difference in means"
##     else "mean"
##     attr(cint, "conf.level") <- conf.level
##     rval <- list(statistic = tstat, parameter = df, p.value = pval, 
##         conf.int = cint, estimate = estimate, null.value = mu, 
##         alternative = alternative, method = method, data.name = dname)
##     class(rval) <- "htest"
##     return(rval)
## }
## <bytecode: 0x412a7c8>
## <environment: namespace:stats>

Tests statistiques

Ce chapitre a pour but de présenter de manière non exhaustive certains tests statistiques employés fréquemment en statistique inférentielle.

Comme on l’a vu précédemment, la détermination des seuils de significativité (p) se fait grâce aux fonctions associées à chaque distribution.

1 - pnorm(167, mean=150, sd=10)
## [1] 0.04456546
1 - pbinom(8,10,0.5)
## [1] 0.01074219

Test du khi-deux

Soit le tableau de contingence A x B suivant à analyser :

A1 A2 A3
B1 13 24 20
B2 10 7 18

Le calcul du test du X^2 associé à ce tableau s’effectue de la manière suivante :

a <- c(13, 24, 20, 10, 7, 18)
chisq.test(matrix(a,2,3,byrow=T))
## 
##  Pearson's Chi-squared test
## 
## data:  matrix(a, 2, 3, byrow = T)
## X-squared = 4.8347, df = 2, p-value = 0.08916

Estimation de la moyenne d’un groupe

L’intervalle de confiance de la moyenne peut être obtenu à l’aide de la fonction t.test :

a <- 10+rnorm(10, sd=10)
t.test(a, conf.level=.01)
## 
##  One Sample t-test
## 
## data:  a
## t = 2.3432, df = 9, p-value = 0.04379
## alternative hypothesis: true mean is not equal to 0
## 1 percent confidence interval:
##  9.327986 9.431148
## sample estimates:
## mean of x 
##  9.379567

Si l’hypothèse de normalité n’est pas soutenable, le test de Wilcoxon (non-paramétrique) peut être utilisé à l’aide de la fonction wilcox.test : ce test des signes permet de déterminer si la médiane du groupe peut être considérée comme significativement différente de 0.

Comparaison de deux groupes

Ce sont les mêmes fonctions - t.test (test paramétrique) et wilcox.test (test non paramétrique) - qui permettent la comparaison entre deux groupes ; dans ce cas, on passe en arguments les deux groupes :

a <- rnorm(10)
b <- rnorm(10, mean=1)
t.test(a, b)
## 
##  Welch Two Sample t-test
## 
## data:  a and b
## t = -1.6056, df = 14.2, p-value = 0.1304
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.0922795  0.2994745
## sample estimates:
##  mean of x  mean of y 
## -0.4645026  0.4318999
wilcox.test(a, b)
## 
##  Wilcoxon rank sum test
## 
## data:  a and b
## W = 25, p-value = 0.06301
## alternative hypothesis: true location shift is not equal to 0
c <- c(a, b)
x <- gl(2, 10, 20)
t.test(c ~ x)
## 
##  Welch Two Sample t-test
## 
## data:  c by x
## t = -1.6056, df = 14.2, p-value = 0.1304
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.0922795  0.2994745
## sample estimates:
## mean in group 1 mean in group 2 
##      -0.4645026       0.4318999
wilcox.test(c ~ x)
## 
##  Wilcoxon rank sum test
## 
## data:  c by x
## W = 25, p-value = 0.06301
## alternative hypothesis: true location shift is not equal to 0

Analyse de variance sur un facteur

Lorsque l’on est en présence d’un ensemble de k observations indépendantes (un seul facteur inter-sujets), on peut comparer leurs moyennes respectives à l’aide de la fonction aov (ou selon un modèle linéaire général, avec la fonction lm).

x<-rnorm(100)
a <- gl(4, 25, 100)
plot(x ~ a)

r <- aov(x ~ a)
anova(r)
## Analysis of Variance Table
## 
## Response: x
##           Df Sum Sq Mean Sq F value Pr(>F)
## a          3  6.175  2.0584  2.0243 0.1156
## Residuals 96 97.617  1.0168
pairwise.t.test(x, a)
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  x and a 
## 
##   1    2    3   
## 2 0.86 -    -   
## 3 0.86 0.11 -   
## 4 0.95 0.95 0.48
## 
## P value adjustment method: holm
t.test(x[a==1], x[a==2])
## 
##  Welch Two Sample t-test
## 
## data:  x[a == 1] and x[a == 2]
## t = 1.2701, df = 47.922, p-value = 0.2102
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1905467  0.8441309
## sample estimates:
##  mean of x  mean of y 
##  0.2104228 -0.1163693

Anova sur deux facteurs

Avec deux facteurs inter-sujets, le principe d’analyse est le même, mais on étudie également l’interaction entre les deux facteurs.

x <- rnorm(100)
a <- gl(2, 50, 100)
b <- gl(2, 25, 100)
plot(x ~ factor(a:b))

interaction.plot(a, b, x)

l <- aov(x~a*b)
anova(l)
## Analysis of Variance Table
## 
## Response: x
##           Df Sum Sq Mean Sq F value Pr(>F)
## a          1  1.630 1.63046  2.3045 0.1323
## b          1  0.000 0.00018  0.0003 0.9874
## a:b        1  1.158 1.15773  1.6364 0.2039
## Residuals 96 67.920 0.70750

Anova sur des protocoles de mesures répétées

Avec un seul facteur intra-sujet, on procèdera ainsi :

subject <- gl(10, 3, 30)
cond <- gl(3, 1, 30)
x <- rnorm(30)
interaction.plot(cond, subject, x)

summary(aov(x ~ cond + Error(subject/cond)))
## 
## Error: subject
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  9  5.504  0.6115               
## 
## Error: subject:cond
##           Df Sum Sq Mean Sq F value Pr(>F)
## cond       2  3.224  1.6122   1.959   0.17
## Residuals 18 14.810  0.8228

Avec deux facteurs intra, la démarche est à peu près identique :

subject <- gl(10,4,40)
cond1 <- gl(2,1,40)
cond2 <- gl(2,2,40)
table(cond1,cond2)
##      cond2
## cond1  1  2
##     1 10 10
##     2 10 10
x <- rnorm(40)
plot(x ~ factor(cond1:cond2))

interaction.plot(cond1,cond2,x)

interaction.plot(cond1,subject,x)

interaction.plot(cond2,subject,x)

summary(aov(x ~ cond1 * cond2 + Error(subject/(cond1 * cond2))))
## 
## Error: subject
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  9  10.66   1.185               
## 
## Error: subject:cond1
##           Df Sum Sq Mean Sq F value Pr(>F)
## cond1      1  4.661   4.661   1.844  0.208
## Residuals  9 22.743   2.527               
## 
## Error: subject:cond2
##           Df Sum Sq Mean Sq F value Pr(>F)
## cond2      1  0.053  0.0527   0.123  0.734
## Residuals  9  3.858  0.4287               
## 
## Error: subject:cond1:cond2
##             Df Sum Sq Mean Sq F value Pr(>F)
## cond1:cond2  1  0.077  0.0772    0.06  0.812
## Residuals    9 11.531  1.2812

Régression linéaire

Comme nous l’avons vu dans le cas des distributions conjointes (cf. section 4.1.2), la démarche pour effectuer de la régression linéaire est la suivante :

a <- rnorm(100)
b <- 2*a + rnorm(100)
plot(b~a)
r <- lm(b ~ a)
anova(r)
## Analysis of Variance Table
## 
## Response: b
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## a          1 460.87  460.87  412.77 < 2.2e-16 ***
## Residuals 98 109.42    1.12                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
abline(r)

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