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

Introduction à l’analyse statistiques des données avec R

Posted on avril 30, 2005 in stats

Introduction à R

R est un logiciel pour l’analyse statistique des données. Il fournit les procédures usuelles (t-tests, anova, tests non paramétriques…) et possède des possibilités graphiques performantes pour explorer les données. Pouvant être utilisé aussi bien en mode interactif qu’en mode batch, R est un logiciel libre, dont le code source est disponible et qui peut être recopié et diffusé gratuitement. Des versions compilées de R sont disponibles pour Linux, Windows et Mac OS X.

Installation du système de base

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

Je recommande vivement d’installer l’environnment 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

Après avoir installé le système de base, vous pouvez installer des modules supplémentaires (‘packages’) qui ajoutent des fonctions à R.

Pour l’analyse des données d’expériences, les paquetages suivants sont intéressants: - 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

Ces modules sont disponibles sur les sites CRAN dans la section Contributed extension packages.

Pour installer un module sous Windows, dans RGui, utiliser le menu ‘Package/Install package from CRAN’ (il faut être connecté à Internet).

Pour installer un module sous Linux, il faut d’abord télécharger le fichier package.tar.gz du CRAN, puis exécuter R CMD INSTALL package.tar.gzen ligne de commande.

On peut également utiliser la fonction install.packages() à l’intérieur de R.

Documentation

De nos jours, beaucoup de gens trouvent naturel de pouvoir utiliser les logiciels sans lire de documentation. Si cela est raisonnable pour les logiciels qui réalisent des opérations assez simples, c’est dangereux avec les logiciels qui effectuent des opérations conceptuellement compliquées. Dans le cas de R, qui comprend de nombreuses commandes, il est illusoire d’envisager utiliser ce logiciel sans lire un minimum de documentation. Notre expérience est que les premières heures d’analyse de données avec R nécessitent de fréquents recours aux documentations, mais lorsqu’on est devenu à l’aise, alors il n’y a pratiquement plus besoin de s’y référer.

Il est donc utile de savoir où chercher l’information a propos de R.

Pour les débutants, on trouve sur Internet un bon nombre de documents sur R, notamment dans la section “Documentation/Contributed” du site www.r-project.org. Mentionnons en particulier :

R possède aussi une documentation officielle, sous forme de fichiers pdf et html, qui est copiée sur votre disque dur lors de l’installation du logiciel. Dans l’interface graphique sous Windows, les manuels au format pdf sont accessibles dans les menus Help/Manuals. Il est fortement conseillé de parcourir, au minimum, les deux documents “An Introduction to R” et “R Data Import et Export”.

Les manuels sont également accessibles sous forme html, dans le menu Help/Html help sous Windows, et en tapant help.start() sous Linux. Cela ouvre votre navigateur Internet sur une page web locale qui contient divers liens, entre autres vers ces manuels. Par exemple, le lien Packages/base liste les commandes de bases de R.

Il existe plusieurs livres publiés qui traitent de R. Pour les débutants, les deux livres suivants peuvent offrir une aide utile :

  • Introductory statistics with R par Peter Dalgaard édité par Springer-Verlag.
  • An R and S-plus companion to applied regression par John Fox, édité par Sage publications.

Pour un niveau plus avancé:

  • Analyzing Linguistic Data par Harald Baayen
  • Regression Modeling Strategies par Franck Harrell
  • Modern Applied Statistics with S-PLUS par Venables et Ripley.
  • Mixed-Effects Models in S and S-PLUS par Pinheiro et Bates

Premiers pas

L’interaction avec R se fait en tapant des commandes dans la fenêtre R Console.

Entrer des commandes dans la console R

Pour commencer, vous pouvez utiliser R comme une calculatrice. Cliquez dans la fenêtre ‘R Console’, puis tapez:

2 + 3
## [1] 5

Poursuivez avec:

    a = 5
    a + 8
## [1] 13

Le principe de R est le suivant : vous entrez une ligne de commande, et quand vous tapez sur ‘Entrée’, R lit cette ligne et effectue l’opération demandée.

Essayez maintenant les commandes suivantes :

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

plot(a, b, pch=16, col=2)

La commande plot provoque l’affichage d’une fenêtre graphique

Cliquez à nouveau dans la fenêtre R Console, puis tapez:

a=c(3, 4, 6, 7, 8, 9)
a
## [1] 3 4 6 7 8 9
length(a)
## [1] 6
b = c('alpha', 'beta')
b
## [1] "alpha" "beta"
length(b)
## [1] 2

La variable ‘a’ contient un vecteur numérique à six éléments.

La variable ‘b’ contient un vecteur contenant deux chaînes de caractères.

Les concepts de vecteur et de variable sont essentiels dans R. On y reviendra plus tard ; pour le moment, retenez que :

  • un vecteur n’est rien d’autre qu’un suite d’items qui ont tous le même type (numérique, chaîne de caractères, …). C’est l’objet de base dans R.
  • une variable contient un objet, et permet de le retrouver sans le ré-écrire en entier.

Comme on l’a déjà vu, la liste des variables 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] 54.82851
sd(a) 
## [1] 12.7692
max(a)
## [1] 76.55416
summary(a)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   27.21   45.36   54.01   54.83   61.84   76.55
hist(a)

boxplot(a)

stripchart(a)

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

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

t.test(x1, x2)
## 
##  Welch Two Sample t-test
## 
## data:  x1 and x2
## t = -5.0761, df = 17.656, p-value = 8.344e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -25.10608 -10.39300
## sample estimates:
## mean of x mean of y 
##  94.37346 112.12300
plot(x1 , x2)

summary(lm(x2 ~ x1))
## 
## Call:
## lm(formula = x2 ~ x1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.3081  -3.9215  -0.4167   5.2390   9.7874 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 110.4242    29.0897   3.796  0.00527 **
## x1            0.0180     0.3072   0.059  0.95471   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.691 on 8 degrees of freedom
## Multiple R-squared:  0.0004291,  Adjusted R-squared:  -0.1245 
## F-statistic: 0.003434 on 1 and 8 DF,  p-value: 0.9547

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")
help.start()

Quitter R

La fenêtre R Console étant active, sélectionnez “File/Exit” et répondez “Oui” à la question “Save workspace image?”

Et voilà…

Tout votre travail est-il perdu ?

Non. 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())”.

Il est possible de choisir le nom de fichier où est sauvegardé le workspace (par défaut “.RData”). Cela permet de faire plusieurs analyses indépendantes sans les mélanger. (Voir les menus File/Load workspace/ Save Workspace). Une alternative plus recommandée et de créer un dossier pour chaque analyse de données indépendantes.

Sauvegarder les commandes dans un script

Tapez la commande history(). Une fenêtre s’affiche listant les dernières commandes que vous avez tapées (voir figure ).

Une manière efficace de travailler avec R consiste à sauvegarder les commandes au fur et à mesure dans un fichier texte. Pour cela, en parallèle avec R, ouvrez un éditeur de fichier texte (le plus simple d’entre eux, bien qu’il soit très limité, est le bloc-notes de Windows disponible dans les accessoires).

En utilisant le copier/coller, copier dans le fichier texte les commandes qui font l’essentiel de l’analyse. A la fin de votre session de travail, sauvez ce fichier avec un nom explicite (par exemple le nom de l’expérience) et une extension “.R”.

Quand vous reprendrez cette analyse quelques jours ou mois plus tard, vous pourrez réutiliser ce fichier, qu’on appelle habituellement un script. R vous permettra de ré-executer les commandes de ce script en utilisant la commande source.

Faites un essai: créez un fichier qui contient les lignes suivantes :

a = rnorm(100)
    b = rnorm(100)
    summary(a)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -2.66900 -0.48600  0.01823  0.10450  0.76300  3.03800
    summary(b)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -2.43000 -0.74310 -0.02984 -0.10030  0.55590  1.99500
    cor.test(a,b)
## 
##  Pearson's product-moment correlation
## 
## data:  a and b
## t = 0.40573, df = 98, p-value = 0.6858
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1567286  0.2354743
## sample estimates:
##        cor 
## 0.04095018

Sauvez-le dans “Mes documents”, sous le nom “test.R”.

Dans R, utilisez le menu “File/Change Dir” pour aller dans “Mes Documents”. Puis tapez:

 source('test.R',echo=T)

Vérifiez que cela marche.

Sous Linux, il n’est pas nécessaire de démarrer R: on peut entrer R BATCH script.R sur une ligne de commande dans un terminal et les résultats sont écrits automatiquement dans le fichier script.Rout.

Sauver les résultats d’une analyse

Les commandes et les résultats des analyses statistiques et les graphiques peuvent être copiés/collés dans un document.

Les résultats (sans les commandes) peuvent être copiés automatiquement dans un fichier texte grâce à “sink”. Tapez:

sink('monanalyse.txt',split=T)
a = 1:10
mean(a)
summary(a)
sink()

Puis ouvrez le fichier “monanalyse.txt”.

Les graphiques peuvent être sauvés directement dans des fichiers graphiques en utilisant les commandes postscript, jpeg ou png (voir l’aide en ligne de ces fonctions).

Mentionnons le paquetage R2HTML qui permet de créer des rapports au format html de façon semi-automatique.

Organisation du travail

L’expérience prouve que la meilleure stratégie est de créer un répertoire (dossier) par analyse de données, et d’y disposer: (a) les fichiers de données brutes; (2) le fichier script contenant les commandes R; (3) le workspace et le(s) fichiers(s) résultats (textes et graphiques).

Manipulations de base

Objets

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

Un type de vecteur particulièrement utile est le type factor. Les facteurs sont des vecteurs utilisés pour classifier les valeurs d’autres vecteurs (les facteurs sont des “variables indicatrices”). Par exemple, étant donné 100 scores provenant de plusieurs groupes de sujets, une variable facteur peut désigner ces sous-groupes.

(a=factor(c(rep('alpha', 10), rep('beta', 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.74619942 -2.03415825 -2.45078493  0.78517866  1.50267863
##  [6] -0.51583617  0.15303998 -0.78983022 -0.36816267  1.96911337
## [11] -0.34784589 -0.39183255  1.27439312 -0.20496299  0.45994364
## [16]  0.04581141  1.15732722  0.42417014  0.90968587  0.05807693
## [21]  1.18935002  1.04772407  0.18628689 -0.76410899 -1.37086541
## [26] -0.20208121  1.48181630 -0.07871632 -0.90276822  0.01289568
## [31]  0.42292397 -0.93048265  0.19481941  0.87724196  1.29420092
## [36]  0.60231022 -1.15562056  0.23188275  0.81197647 -1.35059223
## [41] -0.59152969 -1.99065460 -1.33351722 -0.56438585 -0.48236896
## [46] -0.85680073  0.77597240 -1.08945055
    tapply(x, b, mean)
##          a          b          c 
## -0.2814362 -0.1861379  0.2397781
    boxplot(x ~ b)

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

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

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

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

(a = rnorm(50))
##  [1] -0.49004055 -1.29394884  0.14884944 -0.65218245  0.32989433
##  [6] -2.54890235 -0.16944231 -0.75216359  2.07353704 -1.19347752
## [11] -0.44032492  0.64959852  1.02221616 -1.69837792  1.44322146
## [16]  1.36284979 -0.77778731  0.49949522 -1.15299325 -0.80528923
## [21] -0.81853154 -0.56139744  1.36942334  1.59438794 -1.18980606
## [26]  0.06635126  0.09496587 -0.72871283 -1.60682388  0.58966812
## [31]  0.22358980  1.37768291  1.29782356  0.65707543 -0.13396692
## [36] -1.70621165  0.18842799  0.03155099 -0.04787062 -0.41269713
## [41] -0.65399758 -1.01770287  0.06493552 -0.45009732  2.25663057
## [46]  1.06174241 -0.32123939 -0.23853805 -0.45239721  1.61445166
    a[1]
## [1] -0.4900405
    a[2]
## [1] -1.293949
    a[c(1, 3, 5)]
## [1] -0.4900405  0.1488494  0.3298943
    a>0
##  [1] FALSE FALSE  TRUE FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE
## [12]  TRUE  TRUE FALSE  TRUE  TRUE FALSE  TRUE FALSE FALSE FALSE FALSE
## [23]  TRUE  TRUE FALSE  TRUE  TRUE FALSE FALSE  TRUE  TRUE  TRUE  TRUE
## [34]  TRUE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE  TRUE FALSE
## [45]  TRUE  TRUE FALSE FALSE FALSE  TRUE
    a[a>0]
##  [1] 0.14884944 0.32989433 2.07353704 0.64959852 1.02221616 1.44322146
##  [7] 1.36284979 0.49949522 1.36942334 1.59438794 0.06635126 0.09496587
## [13] 0.58966812 0.22358980 1.37768291 1.29782356 0.65707543 0.18842799
## [19] 0.03155099 0.06493552 2.25663057 1.06174241 1.61445166
    (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] -0.4900405 -1.2939488  0.1488494 -0.6521825  0.3298943 -2.5489023
##  [7] -0.1694423 -0.7521636  2.0735370 -1.1934775 -0.4403249  0.6495985
## [13]  1.0222162 -1.6983779  1.4432215  1.3628498 -0.7777873  0.4994952
## [19] -1.1529932 -0.8052892 -0.8185315 -0.5613974  1.3694233  1.5943879
## [25] -1.1898061

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 s’avère très utile pour 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’.

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

Arrays, listes et data.frames

D’autres objets de R sont les les listes, les arrays (vecteurs multidimensionnels) et les data.frames.

Les data.frames sont des listes de vecteurs qui ont tous la même longueur. Les data.frames sont très bien adaptés pour stocker des données présentées sous forme de tableau bi-dimensionnel.

(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.99414613
## 2  1  2 -2.07587886
## 3  1  3 -0.59902787
## 4  1  4  0.57858822
## 5  1  5 -1.75878559
## 6  2  6  0.02992957
## 7  2  7 -1.53670296
## 8  2  8 -1.09951551
## 9  2  9 -0.26293485
## 10 2 10  2.04475841
    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.99414613 -2.07587886 -0.59902787  0.57858822 -1.75878559
##  [6]  0.02992957 -1.53670296 -1.09951551 -0.26293485  2.04475841
    c[1:2,]
##   a b          x
## 1 1 1 -0.9941461
## 2 1 2 -2.0758789

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"    "freq" "mots" "v"    "x"    "x1"   "x2"
    rm(a)
    ls()
## [1] "b"    "c"    "freq" "mots" "v"    "x"    "x1"   "x2"

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.18172575 -0.17663704  0.03035259  0.27473636  1.24954705
##  [6] -1.48911712 -0.36496406 -1.66550503  0.78051965  0.59428103
    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.74619942 -2.03415825 -2.45078493  0.78517866  1.50267863
##  [6] -0.51583617  0.15303998 -0.78983022 -0.36816267  1.96911337
## [11] -0.34784589 -0.39183255  1.27439312 -0.20496299  0.45994364
## [16]  0.04581141  1.15732722  0.42417014  0.90968587  0.05807693
## [21]  1.18935002  1.04772407  0.18628689 -0.76410899 -1.37086541
## [26] -0.20208121  1.48181630 -0.07871632 -0.90276822  0.01289568
## [31]  0.42292397 -0.93048265  0.19481941  0.87724196  1.29420092
## [36]  0.60231022 -1.15562056  0.23188275  0.81197647 -1.35059223
## [41] -0.59152969 -1.99065460 -1.33351722 -0.56438585 -0.48236896
## [46] -0.85680073  0.77597240 -1.08945055
    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 et read.table 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.086277833 0.004575344 0.923430190 0.752197812 0.553001953
##  [6] 0.797581331 0.448278308 0.238572420 0.162277250 0.464658049
    rnorm(10) # distribution normale
##  [1] -1.8401052 -0.2653521  2.9060628  0.2331750 -0.1691501 -0.2623380
##  [7] -1.3281432  1.4255977 -0.9852099 -1.1241900
    rnorm(10,mean=100)
##  [1] 100.99702 100.99634 100.65538  99.76738 101.17303 101.96837 100.79886
##  [8] 101.22319 100.14854 100.36403
    rbinom(10,size=1,prob=.5) # distribution binomiale
##  [1] 0 0 1 1 0 1 0 1 1 1

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.4292  0.9306  2.0340  2.1440  2.9900  4.8400
    s2=rnorm(100,mean=2)
    summary(s2)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -1.142   1.430   2.049   2.100   2.710   5.691
    s3=rnorm(10000,mean=2)
    summary(s3)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -1.945   1.332   2.011   1.999   2.670   5.897
    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)
## [1] 0.9772499
    pt(3,df=10)             # fonction de répartition de la loi du t de Student
## [1] 0.9933282
    qnorm(.99)              # donne la valeur associée au 99ème centile d'une distribution normale
## [1] 2.326348
    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.08516488
    sd(a) # écart-type corrigé
## [1] 1.007532
    summary(a)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -3.53300 -0.63330 -0.09532 -0.08516  0.49700  2.33200
    boxplot(a)

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

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.001136598
    clmean(a)
## [1] -0.07499423

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: 0x2f10d78>
## <environment: namespace:base>
    t.test
## function (x, ...) 
## UseMethod("t.test")
## <bytecode: 0x376e860>
## <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: 0x4e8f7d0>
## <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 = 1.9548, df = 9, p-value = 0.08233
## alternative hypothesis: true mean is not equal to 0
## 1 percent confidence interval:
##  6.923750 7.015638
## sample estimates:
## mean of x 
##  6.969694

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 = -4.3676, df = 16.671, p-value = 0.0004375
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.426473 -0.844170
## sample estimates:
##  mean of x  mean of y 
## -0.3635432  1.2717781
    wilcox.test(a,b)
## 
##  Wilcoxon rank sum test
## 
## data:  a and b
## W = 8, p-value = 0.0007253
## 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 = -4.3676, df = 16.671, p-value = 0.0004375
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.426473 -0.844170
## sample estimates:
## mean in group 1 mean in group 2 
##      -0.3635432       1.2717781
    wilcox.test(c ~ x)
## 
##  Wilcoxon rank sum test
## 
## data:  c by x
## W = 8, p-value = 0.0007253
## 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  0.707 0.23561  0.2621 0.8525
## Residuals 96 86.285 0.89880
    pairwise.t.test(x, a)
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  x and a 
## 
##   1 2 3
## 2 1 - -
## 3 1 1 -
## 4 1 1 1
## 
## 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 = 0.13762, df = 44.712, p-value = 0.8912
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.4850950  0.5562323
## sample estimates:
##   mean of x   mean of y 
## -0.09051978 -0.12608843

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   0.682 0.68177  0.5961 0.4420
## b          1   0.116 0.11622  0.1016 0.7506
## a:b        1   2.563 2.56262  2.2406 0.1377
## Residuals 96 109.796 1.14371

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  10.03   1.115               
## 
## Error: subject:cond
##           Df Sum Sq Mean Sq F value Pr(>F)
## cond       2  1.997  0.9985   0.581   0.57
## Residuals 18 30.956  1.7198

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  13.64   1.516               
## 
## Error: subject:cond1
##           Df Sum Sq Mean Sq F value Pr(>F)  
## cond1      1  3.250    3.25   3.532 0.0929 .
## Residuals  9  8.282    0.92                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: subject:cond2
##           Df Sum Sq Mean Sq F value Pr(>F)
## cond2      1  0.496  0.4957   0.365  0.561
## Residuals  9 12.215  1.3572               
## 
## Error: subject:cond1:cond2
##             Df Sum Sq Mean Sq F value Pr(>F)  
## cond1:cond2  1  5.586   5.586   5.243 0.0478 *
## Residuals    9  9.590   1.066                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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 469.54  469.54  468.27 < 2.2e-16 ***
## Residuals 98  98.26    1.00                      
## ---
## 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))

    summary(lm(c~a*b))
## 
## Call:
## lm(formula = c ~ a * b)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.21832 -0.64111  0.04496  0.68300  2.33801 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.159704   0.120964   1.320    0.190    
## a            5.052426   0.248680  20.317   <2e-16 ***
## b           -0.103677   0.109639  -0.946    0.347    
## a:b         -0.003639   0.037842  -0.096    0.924    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.037 on 96 degrees of freedom
## Multiple R-squared:  0.9527, Adjusted R-squared:  0.9513 
## F-statistic: 645.2 on 3 and 96 DF,  p-value: < 2.2e-16

Exemples d’analyses de données

Ces exemples proviennent principalement du site web “Analyse Statistique des Données en Psychologie (ASDP)” de l’UFR de Psychologie de l’université Paris 5 (piaget.psycho.univ-paris5.fr/, lien “Analyse des Données” puis “Données”).

Dossier sommeil

Lors d’une expérimentation médicale, on a relevé le temps de sommeil T de 10 patients (facteur Sujet) sous l’effet de deux médicaments (d’où le facteur Médicament M). Chaque sujet a pris successivement l’un et l’autre des deux médicaments.

Source  

Student (1908) The probable error of a mean, Biometrika, VI, 1-25.

Données  

Fichier sommeil.txt

Question  

Ces données ont été recueillies pour tester l’hypothèse que le médicament m2 est plus efficace que le médicament m1. Est-ce le cas ?

Une solution   Voir l’exemple de script listé en A.1 et B.1 (Statistica).

Dossier pédagogie

Lors d’une expérimentation pédagogique, on désire comparer l’efficacité de quatre méthodes d’enseignements.

On dispose des notes obtenues à un examen par quatre groupes d’élèves ayant chacun reçu un des 4 types d’enseignements.

Source:  

Données fictives.

Données  

Fichier pedago.txt

Questions  

Comparer les résultats obtenus en fonction des méthodes.

Une solution   Voir l’exemple de script listé en A.2 et B.2 (Statistica).

Dossier négligence

Une recherche a porté sur la “pseudo-négligence” qu’on observe chez des sujets normaux. Ce nom provient des similarités qu’elle présente avec l’hémi-négligence (atteinte de la moitié du champ visuel) de sujets atteints d’une lésion cérébrale. La tâche des sujets consiste à déterminer le milieu subjectif d’une baguette de 24cm avec la seule aide d’informations kinesthésiques. La pseudo-négligence se traduit par une déviation systématique vers la droite (pour les droitiers) de ce milieu subjectif par rapport au milieu objectif de la baguette.

Les données portent sur 24 femmes droitières (facteur S) réparties selon 2 conditions (12 sujets pour chacune): active (c1) où le sujet peut librement déplacer son doigt posé sur un curseur mobile le long de la baguette; ou passive (c2) où le sujet commande un moteur déclenchant le mouvement de la baguette dans un sens ou dans l’autre, alors que son doigt ne bouge pas (facteur C). Chaque sujet exécute cette tâche dans 6 situations expérimentales obtenues par le croisement de: la main utilisée, gauche (m1) ou droite (m2); et l’orientation du regard, 30° à gauche (o1), 0° (o2) ou 30° à droite (o3) (facteurs M et O). Pour chaque sujet et chaque situation on mesure la déviation en cm entre le milieu subjectif et le milieu objectif de la baguette. Une déviation à droite est notée par une valeur positive, à gauche par une valeur négative.

On s’intéresse ici à l’effet de la condition (C) lorsque le sujet utilise sa main habituelle (m2) (Rappel : tous les sujets sont droitiers) et lorsqu’il se trouve en face du milieu de la baguette (avec l’orientation à 0 degrés)

Source  

Chokron, Imbert (1993) - Egocentric reference and asymmetric perception of space, Neuropsychologia, 31, 3, 267-275. D’après J.M. Bernard (1994) - Structure des données, d

Données   fichier neglige2.txt

Questions  

Importer ces données, les visualiser, comparer les groupes. Conclusion?

Une solution  

Voir l’exemple de script listé en A.3 et B.3 (Statistica).

Dossier family

Étude réalisée au USA sur les origines des stéréotypes liés au sexe. 35 familles choisies au hasard et ayant une fille ainée (ou fille unique) en “ninth grade” (Troisième).

Le père a répondu à un questionnaire sur ses intérêts pour le sport, noté sur une échelle numérique de 0 à 50 (FATH)

La mère a répondu au même questionnaire (MOTH)

Le professeur d’éducation physique de chacune des filles a noté les performances physiques générales de la fille de 0 à 20 (PROF).

La fille a répondu également au questionnaire d’inérêt pour le sport (GIRL).

Source   Hays, W.L. (1994) - Statistics, Fort Worth: Harcourt Brace College Publishers (5ème édition), p.671-672

Données   family.txt

Questions   Que faire avec ces données?

Une solution   Voir l’exemple de script listé en A.4 et B.4 (Statistica).

Dossier IO

En 1980, on a interrogé des lycéens (garçons et filles) sur leurs intentions d’orientation après le bac (études scientifiques, littéraires ou techniques).

Source   Il s’agit de données en partie fictives, inspirées d’un exemple de M. Reuchlin.

Données   Fichier io.txt

Questions   Peut-on dire que l’orientation envisagée est liée au sexe chez l’ensemble des lycéens de cette année 1980 ?

Une solution   Voir l’exemple de script listé en A.5 et B.5 (Statistica).

Solutions sous R

Nous proposons ici des scripts pour analyser les exemples du chapitre 6. Il y a plusieurs manières de résoudre le même problème avec R. Par conséquent, vos scripts peuvent différer.

Dossier sommeil

sommeil<-read.table('donnees/sommeil1.txt', header=TRUE)
    sommeil
##      M1   M2
## S1  5.7  6.9
## S2  3.4  5.8
## S3  4.8  6.1
## S4  3.8  5.1
## S5  4.9  4.9
## S6  8.4  9.4
## S7  8.7 10.5
## S8  5.8  6.6
## S9  5.0  9.6
## S10 7.0  8.4
    attach(sommeil)                                                                                  
    summary(M1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   3.400   4.825   5.350   5.750   6.700   8.700
    summary(M2)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.900   5.875   6.750   7.330   9.150  10.500
    plot(M1, M2, xlim=c(0,10), ylim=c(0,10), col=2)
    identify(M1, M2, sommeil)
## integer(0)
    abline(0, 1)

    stripchart(M2 - M1, method='stack')

    t.test(M2 - M1)
## 
##  One Sample t-test
## 
## data:  M2 - M1
## t = 4.0621, df = 9, p-value = 0.002833
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.7001142 2.4598858
## sample estimates:
## mean of x 
##      1.58
    t.test(M1, M2, paired=T)
## 
##  Paired t-test
## 
## data:  M1 and M2
## t = -4.0621, df = 9, p-value = 0.002833
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.4598858 -0.7001142
## sample estimates:
## mean of the differences 
##                   -1.58
    detach()

Dossier pedago

a<-read.table('donnees/pedago.txt', header=TRUE)
    attach(a)
    boxplot(notes ~ pedago)

    stripchart(notes~pedago,method='stack',vertical=T)

    tapply(notes,pedago,mean)
##        a        b        c        d 
## 14.88000 14.00000 16.53571 16.27778
    tapply(notes,pedago,sd)
##        a        b        c        d 
## 2.088061 1.777047 2.151350 2.739806
    tapply(notes,pedago,summary)
## $a
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   10.00   14.00   15.00   14.88   17.00   19.00 
## 
## $b
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      11      13      14      14      15      17 
## 
## $c
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   11.00   15.00   16.50   16.54   18.00   20.00 
## 
## $d
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   13.00   14.25   16.00   16.28   17.50   22.00
    barplot(t(tapply(notes,pedago,mean)))

    m<-aov(notes~pedago)
    summary(m)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## pedago       3   95.5   31.82   6.636 0.000434 ***
## Residuals   87  417.2    4.80                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    TukeyHSD(m)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = notes ~ pedago)
## 
## $pedago
##           diff        lwr      upr     p adj
## b-a -0.8800000 -2.6008450 0.840845 0.5405529
## c-a  1.6557143  0.0773421 3.234086 0.0360261
## d-a  1.3977778 -0.3753845 3.170940 0.1728758
## c-b  2.5357143  0.8563415 4.215087 0.0008819
## d-b  2.2777778  0.4141419 4.141414 0.0101278
## d-c -0.2579365 -1.9908790 1.475006 0.9797521
    plot(TukeyHSD(m))

Dossier négligence

d<-read.table('donnees/neglige4.txt')
    x<-d$V1
    a<-gl(2,12,24)
    b<-gl(2,6,24)
    table(a,b)
##    b
## a   1 2
##   1 6 6
##   2 6 6
    tapply(x,list(a=a,b=b),mean)
##    b
## a       1         2
##   1 -0.25 -1.200000
##   2 -0.75  1.283333
    interaction.plot(a,b,x)

    l<-aov(x~a*b)
    summary(l)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## a            1   5.90   5.900   9.557 0.005756 ** 
## b            1   1.76   1.760   2.851 0.106834    
## a:b          1  13.35  13.350  21.623 0.000154 ***
## Residuals   20  12.35   0.617                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    model.tables(l,se=T)
## Tables of effects
## 
##  a 
## a
##       1       2 
## -0.4958  0.4958 
## 
##  b 
## b
##        1        2 
## -0.27083  0.27083 
## 
##  a:b 
##    b
## a   1       2      
##   1  0.7458 -0.7458
##   2 -0.7458  0.7458
## 
## Standard errors of effects
##              a      b    a:b
##         0.2268 0.2268 0.3208
## replic.     12     12      6
    t.test(x[a==1 & b==1],x[a==1 & b==2])
## 
##  Welch Two Sample t-test
## 
## data:  x[a == 1 & b == 1] and x[a == 1 & b == 2]
## t = 3.3075, df = 7.7985, p-value = 0.01113
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.2846604 1.6153396
## sample estimates:
## mean of x mean of y 
##     -0.25     -1.20
    t.test(x[a==2 & b==1],x[a==2 & b==2])
## 
##  Welch Two Sample t-test
## 
## data:  x[a == 2 & b == 1] and x[a == 2 & b == 2]
## t = -3.5444, df = 9.9716, p-value = 0.005341
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.3120698 -0.7545969
## sample estimates:
## mean of x mean of y 
## -0.750000  1.283333

Dossier family

fam<-read.table('donnees/family.txt', header=T)
    fam
##     FATH MOTH PROF GIRL
## I01   25    8   24   30
## I02   32   30   13   30
## I03   25   25   15   25
## I04   36   31   10   35
## I05   37   36    9   39
## I06   31   33   10   30
## I07   28   22    8   27
## I08   26   18    5   14
## I09   20   14   13   39
## I10   35   30   15   42
## I11   42   36   11   41
## I12   40   30   11   38
## I13   32   27   12   39
## I14   33   38   10   37
## I15   45   40   18   47
## I16   30   31   13   30
## I17   29   25   12   36
## I18   28   24    4   12
## I19   24   27   10   26
## I20   18   20    8   20
## I21   37   35   16   43
## I22   27   29   15   31
## I23   39   32   14   40
## I24   34   30    7   31
## I25   48   40   12   36
## I26   37   31    6   21
## I27   44   41   19   44
## I28   23   26    5   11
## I29   29   24    8   27
## I30   38   30    9   24
## I31   41   20   15   40
## I32   39   31    7   32
## I33   27   29    6   10
## I34   48   43   14   37
## I35   33   30    7   19
    attach(fam)
    data<-as.matrix(fam[,-1])
    pairs(data,panel=panel.smooth)

    cor(data)
##            MOTH       PROF      GIRL
## MOTH 1.00000000 0.03604651 0.4084647
## PROF 0.03604651 1.00000000 0.6732057
## GIRL 0.40846468 0.67320575 1.0000000
    cor.test(FATH, GIRL)
## 
##  Pearson's product-moment correlation
## 
## data:  FATH and GIRL
## t = 4.2286, df = 33, p-value = 0.000175
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3234750 0.7732961
## sample estimates:
##       cor 
## 0.5928176
    cor.test(MOTH, GIRL)
## 
##  Pearson's product-moment correlation
## 
## data:  MOTH and GIRL
## t = 2.5707, df = 33, p-value = 0.01485
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.08707006 0.65284617
## sample estimates:
##       cor 
## 0.4084647
    cor.test(PROF, GIRL)
## 
##  Pearson's product-moment correlation
## 
## data:  PROF and GIRL
## t = 5.2299, df = 33, p-value = 9.363e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4382858 0.8220343
## sample estimates:
##       cor 
## 0.6732057
    l<-lm(GIRL ~ FATH + MOTH + PROF)
    summary(l)
## 
## Call:
## lm(formula = GIRL ~ FATH + MOTH + PROF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.390  -4.947   1.108   3.962  14.095 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -3.3967     5.1153  -0.664   0.5116    
## FATH          0.4605     0.2229   2.066   0.0473 *  
## MOTH          0.1571     0.2184   0.720   0.4771    
## PROF          1.2994     0.2547   5.103  1.6e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.141 on 31 degrees of freedom
## Multiple R-squared:  0.6493, Adjusted R-squared:  0.6153 
## F-statistic: 19.13 on 3 and 31 DF,  p-value: 3.295e-07
    detach(fam)

Dossier IO

a <- read.table('donnees/io90.txt', header=T)
chisq.test(a)
## 
##  Pearson's Chi-squared test
## 
## data:  a
## X-squared = 6.6667, df = 2, p-value = 0.03567