Cours 7
Les régressions linéaires
Les analyses statistiques
avec R
Modèles linéaires ,lm()
modèles linéaires généralisés,glm()
analyse de variance, aov()
Généralités
• L ’argument principal est une formule du
type:
réponse ~ prédicteur
exemples: data(I=InsectSprays)
aov(sqrt(count) ~ spray,data=I)
équivalent à :
aov(sqrt(I$count) ~ I$spray)
ou à aov(sqrt(I[,1]) ~ I[,2])
Les formules
• y~model ou y est la réponse analysée et model
est un ensemble de termes pour lesquels les
paramètres sont estimés
• attention: les symboles arithmétiques ont ici
une signification particulière
• exemples:
y~x1+x2 désigne le modèle y=ax1+bx2+c
y~I(x1+x2) désigne le modèle y=a(x1+x2)+c
y~poly(x,2) désigne le modèle y=ax^2+bx+c
y~x1+x2 désigne le modèle y=ax1+bx2+c
y~x1-1 désigne le modèle y=ax1
Les fonctions génériques
• Les objets qui contiennent les résultats
d ’une analyse, ont un attribut particulier,
la classe.
• Certaines fonctions, dites génériques,
permettent d ’extraire des informations d ’un
objet résultat
• exemples: summary()qui a une action différente
sur un objet de classe lm() , aov(),...
apropos("^summary")
[1] "[Link]" "[Link]"
"[Link]"
[4] "[Link]" "[Link]"
"[Link]"
[7] "[Link]" "[Link]"
"summaryRprof"
[10] "summary" "[Link]"
"[Link]"
[13] "[Link]" "[Link]"
"[Link]"
[16] "[Link]" "[Link]"
"[Link]"
[19] "[Link]"
Les fonctions génériques
• plot()
• print() résumé succint
• summary() résumé détaillé
• [Link] nbre de ddl résiduels
• coef coefficients estimés
• residuals résidus
• deviance déviance
• fitted valeurs ajustées par le modèle
• logLik logarithme de la vraisemblance...
Un objet-résultat, produit par
aov(), lm()….
est généralement une liste,bien qu ’il ne soit
pas affiché comme tel,
dont les éléments peuvent être affiché par la
fonction names()
ex names([Link])pour une régréssion linéaire
"coefficients" "residuals" "effects"
"rank" "[Link]" "assign"
"qr" "[Link]" "xlevels" "call"
"terms" "model"
Analyses supplémentaires
à partir d ’un objet
add1():teste successivement tous les termes qui
peuvent être ajoutés à un modèle
drop1():teste successivement tous les termes qui
peuvent être enlevés à un modèle
anova(): calcule une table d ’analyse de
variance ou de deviance pour un ou plusieurs
modèles
predict(): calcule les valeurs prédites pour des
nouvelles données
update(): réajuste un modèle ...
Régréssion linéaire simple
• Sur des données fictives:
x=1:100
x=sample(x,30,replace=TRUE)
x
[1] 62 18 9 67 43 38 57 12 41 29 69 76 77 46 42 75 32 74 6 40 51 88 61 3 38
[26] 71 81 76 94 34
y=3+7*x+rnorm(30,0,100)
Y
[1] 340.61710 254.86969 54.52298 463.78335 379.30676 177.27873 555.98297
[8] -13.48922 273.11081 187.46739 439.59869 380.92303 537.40362 414.12641
[15] 299.09269 494.05965 415.9
plot(x,y)
[Link]=lm(y~x);
Call:
• lm(formula = y ~ x)
Coefficients:
(Intercept) x
13.22 6.71
Droite de régression:
y=6,71 *x +13.22
plot(x,y);
abline([Link])
summary([Link])
• Residuals:
Min 1Q Median 3Q Max
-142.29 -58.64 -17.17 63.33 187.99
• Coefficients:
• Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.2213 37.0548 0.357 0.724
x 6.7105 0.6587 10.188 6.37e-11
• Residual standard error: 90.65 on 28 degrees of
freedom
• Multiple R-Squared: 0.7875, Adjusted R-squared:
0.78
Plus sophistiquée...
• La croissance d ’une bactérie (par jour),
modélisé par N=N0 e^kt
t=2:12;
N=c(55,90,135,245,403,665,1100,1810,3000,4450,
7350);
Le modèle est le suivant;
• t=c(2:12);N=c(55,90,135,245,403,665,1100
,1810,3000,4450,7350)
• T=[Link](t,N,y=log(N));T;
> T
t N y
1 2 55 4.007333
2 3 90 4.499810
3 4 135 4.905275
4 5 245 5.501258…..
Calcul de moyenne et
écart-type
• apply(T,2,mean);
t N y
7.000000 1754.818182 6.475094
• apply(T,2,sd);
t N y
3.316625 2326.625317 1.640357
plot(T$t,T$N)
plot(T$t,T$y)
droite de regression
• ll=lm(y~t,data=T);ll;
Call:
lm(formula = y ~ t, data = T)
Coefficients:
(Intercept) t
3.0142 0.4944
abline(ll);
summary(ll)
• Call:
lm(formula = y ~ t, data = T)
• Residuals:
Min 1Q Median 3Q Max
-0.08656 -0.02117 0.01500 0.02912 0.04802
• Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.014162 0.032947 91.49 1.13e-14 ***
t 0.494419 0.004289 115.27 1.41e-15 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
summary(ll) suite
Residual standard error: 0.04499 on 9 degrees of freedom
Multiple R-Squared: 0.9993, Adjusted R-squared: 0.9992
F-statistic: 1.329e+04 on 1 and 9 DF, p-value: 1.413e-15
Régression linéaire
multiple
• Exemples:
Les tests statistiques
les test du chi-deux
La fonction
[Link](x,y,logical,p)
• Premier exemple:
on lance un dé 300 fois et on obtient le résultat suivant:
1 2 3 4 5 6
43 49 56 45 66 41
x=c(43, 49, 56, 45, 66, 41)
prob=rep(1/6,6)
[Link](x,p=prob)
● Chi-squared test for given probabilities
• data: x
• X-squared = 8.96, df = 5, p-value = 0.1107
Second exemple sur un tableau de
contingence
Exemple d ’un tableau donnant la cécité en fonction
du sexe:
tab=matrix(c(442,514,38,6),nrow=2,byrow=TRUE)
colnames(tab)=c("homme","femme")
rownames(tab)=c("voyant","aveugle")
homme femme
voyant 442 514
aveugle 38 6
X2=[Link](tab,correct=FALSE)
On teste s ’il y a une relation entre
sexe et cécité (l ’hypothèse par défaut
est celle d ’indépendance)
Pearson's Chi-squared test
data: tab
X-squared = 27.1387, df = 1, p-value = 1.894e-07
attributes(x2)
$names
[1] "statistic" "parameter" "[Link]"
"method" "[Link]" "observed"
[7] "expected" "residuals"
$class
[1] "htest »
par exemple:
x2$expected
homme femme
voyant 458.88 497.12
aveugle 21.12 22.88
valeurs attendues sous hypothèse d ’indépendance
x2$residuals
homme femme
voyant -0.787994 0.7570801
aveugle 3.673039 -3.5289413
sum(x2$residuals^2)
27.13874 la somme des carrés des résidus est la
valeur du chi-deux
• Soit le tableau de contingence suivant:
• roux blond brun
• bleu 13 20 7
• marron 24 10 18
• le test du chi-deux d ’indépendance
s ’effectue ainsi:
• [Link](m)
• Pearson's Chi-squared test
• data: m
• X-squared = 10.0494, df = 2, p-value =
0.006574
on teste l ’hypothèse nulle suivante
« H0:il y a indépendance entre la couleur des
yeux et celle des cheveux »
Test sur une moyenne:
[Link]()
• Pour comparer deux ou plusieurs
proportions sur des échantillons de
grande taille:
[Link](x, n, p = NULL,alternative =
c("[Link]", "less", "greater"),
[Link] = 0.95, correct = TRUE)
Classification
La fonction kmeans()
Méthodes de type « nuées dynamiques »:
on suppose que les individus sont des points de R^n
muni de la distance euclidienne
On recherche un partitionnement des individus en classes
• Variance des classes minimale (classes homogènes)
• Variance entre classes maximale
• Rappel : théorème de Huygens: la somme de l ’inertie
interclasse et de l ’inertie intraclasse est constante
(inertie = moyenne des carrés des distances au centre de
gravité)
suite
• Exemple en dimension 2 sur le dataframe D
PROB LECT CARR OPER SYNO PFB SUITE ANAL D70 T M s
1 4.52 1.92 0.93 1.94 1.20 0.08 2.05 1.99 1.37 16.01 1.78 1.21
2 1.50 0.77 0.43 -0.32 1.20 1.26 2.05 1.99 2.30 11.19 1.24 0.85
3 1.77 0.19 1.94 1.19 1.20 0.59 2.05 1.50 1.37 11.80 1.31 0.61
4 1.22 0.19 2.19 0.43 1.20 0.92 2.05 0.51 0.75 9.10 1.01 0.76
5 1.22 1.92 -0.83 1.56 0.64 0.08 1.76 1.01 0.44 7.82 0.87 0.88
6 0.68 0.39 -0.83 -0.32 1.20 -1.77 1.76 0.76 0.18 1.70 0.19 1.08
7 2.05 0.77 0.68 0.43 0.36 1.60 1.76 0.51 1.06 9.23 1.03 0.63
8 0.42 0.39 -0.57 1.56 0.36 -0.93 1.76 0.26 0.44 2.86 0.32 0.90
9 -0.15 -0.19 1.44 -0.32 0.64 0.76 1.76 0.02 1.37 5.34 0.59 0.79
10 -0.42 -0.19 -0.57 1.19 0.64
E=D[,c(1,2)]
cl <- kmeans(E, 4, 20) (donne 4 sous-nuages)
plot(E, col = cl$cluster) (tracé pour un objet de type résultat
de la fct kmeans)
points(cl$centers, col = 1:4, pch = 8)
Autre exemple
• data(airquality); a=airquality [,3:4]
• cl=kmeans(a,3,20)
• plot(a, col = cl$cluster) ;points(cl$centers, col = 1:4, pch
= 8)
le résultat de la fct kmeans est une liste
contenant :les composants suivants:
• cluster: un vecteur d’entiers indiquant la
partition ( le sous-nuage) à laquelle est
affecté chaque point
• centers: la matrice des centres des sous-
nuages
• withinss: The within-cluster sum of squares
for each cluster.
• size: le nombre de points dans chaque sous-
nuage
La fonction hclust()
●
La fonction cutree()
●
