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),

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…).

Manipulation des distributions de probabilités

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

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)

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)

?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')

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))

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

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

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

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

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>

Test du khi-deux

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

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))

summary(lm(c ~a * b))

## 
## Call:
## lm(formula = c ~ a * b)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.33407 -0.58125  0.00417  0.66091  2.12656 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.09815    0.11412  -0.860    0.392    
## a            4.99079    0.18874  26.442   <2e-16 ***
## b            0.05983    0.09138   0.655    0.514    
## a:b          0.02373    0.03095   0.767    0.445    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9213 on 96 degrees of freedom
## Multiple R-squared:  0.974,  Adjusted R-squared:  0.9732 
## F-statistic:  1198 on 3 and 96 DF,  p-value: < 2.2e-16

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).

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