TP1 Corr

Télécharger au format pdf ou txt
Télécharger au format pdf ou txt
Vous êtes sur la page 1sur 6

Sam Meyer - 4BB-BIOSTAT2 2020-2021

TP1 : Bases de la régression linéaire

1 Un exemple complet d’application de la régression


1.1 Construction du modèle
> c<-read.table("conduct.txt")
> par(mfrow=c(1,2)) # réaliser deux graphes cote a cote
> plot(c$V1, c$V2)
> l=lm(c$V2~c$V1)
> abline(l$coefficients)
> plot(c$V1,l$residuals)
> abline(h=0,lt="dashed")

● ●
10

● ●
2



● ● ●
● ● ●

● ● ●
8

● ● ● ●
1


● ●
● ●
● ●
l$residuals

● ● ● ● ●
● ● ●
● ● ●
6
c$V2

● ●

● ●
0

● ●
● ● ●
● ● ● ●
● ●
● ● ● ●
4

● ● ●
● ● ● ●
● ● ● ●
● ● ● ● ● ●
−1

● ● ●

2

● ● ● ● ● ●
● ●

● ●
−2

● ●
0

2 4 6 8 10 2 4 6 8 10

c$V1 c$V1

1. D’après le cours de chimie, pour une solution simple, la conductivité est proportionnelle à la concen-
tration du soluté. Mais en cas de solution complexe (plusieurs solutés), la dépendance peut être plus
complexe et donner lieu à des interactions. Enfin, la mesure peut être affectée par des erreurs d’expéri-
mentation, p. ex. si le zéro n’a pas été effectué correctement et qu’on a un “bruit de fond”.
2. On constate bien une dépendance apparemment linéaire.
3. On suppose y = β0 + βx +  avec  = N (0, σ). Ici, β0 et β sont les deux paramètres de la partie
explicative, et σ est l’écart-type des variations dues aux facteurs non contrôlés.

1.2 Analyse des résultats


6. Il y a autant de résidus que de points, donc ici 50. Les deux paramètres explicatifs sont l’ordonnée à
l’origine et la pente, de valeurs respectives suivantes :
> summary(l)

Call:
lm(formula = c$V2 ~ c$V1)

Residuals:
Min 1Q Median 3Q Max

1
-2.0166 -0.8328 -0.2537 0.8387 2.1041

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.4828 0.3146 1.535 0.131
c$V1 0.9150 0.0507 18.048 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.03 on 48 degrees of freedom


Multiple R-squared: 0.8716, Adjusted R-squared: 0.8689
F-statistic: 325.7 on 1 and 48 DF, p-value: < 2.2e-16

> l$coefficients

(Intercept) c$V1
0.4827680 0.9149489

7. Voir code
8. On voit que la dispersion verticale des points est de l’ordre de 2-3, donc ça indique une valeur de
σ de l’ordre de l’unité. R nous fournit une estimation de 1.03 (residual standard error) calculée sur
50 − 2 = 48 degrés de liberté, donc ça correspond.
9. On voit que les groupes de résidus semblent centrés en 0 et de même variance approximative. Par contre
il est plus difficile de vérifier la normalité sur ce graphe.
10. Oui, puisque l’hypothèse nulle H0 : β = 0 est testée par R (sous réserve de la normalité des résidus) à
l’aide d’un test de Student, qui donne une p-value inférieure à 2.10−16 , ce qui est hautement significatif
quelle que soit la valeur de α raisonnablement choisie. Cela confirme l’impression évidente que nous
avons sur la première figure, puisque cette hypothèse nulle consiste à dire que la pente apparente du
graphe est dûe à un placement aléatoire des points, ce qui est invraisemblable.

1.3 Validation graphique du modèle


Observez les quatre graphes dans R en tapant la commande plot(l) puis Enter.

1.4 Erreurs sur les paramètres et modèle alternatif


11. Le cours nous dit que le nombre β (pente réelle) peut être estimé à l’aide d’une loi de Student à 48
degrés de liberté, en sachant la moyenne b = 0.9150 et son écart-type expérimental 0.0507. Le nombre
β−0.9150
0.0507 a donc 95% de chances d’être situé dans la partie centrale de la loi de Student, entre les bornes
données par les quantiles qt(0.025,48) et qt(0.975,48), qui valent approximativement −2.01 et
2.01. Vous avez le droit de faire l’approximation ±2 pour une loi de Student à partir de 5-10 degrés
de liberté, à condition de bien préciser la loi exacte. En isolant le nombre β, on obtient β ∈ 0.9150 ±
2.01 × 0.0507 = 0.92 ± 0.10 en ne gardant que deux chiffres significatifs vu la largeur de l’intervalle
de confiance.

12. D’après le cours, la droite passe par ȳ au point x̄, et l’écart-type est alors donné par s/ n. En prenant
toujours les quantiles de la loi de Student, on a l’intervalle :
> x=c$V1; xm=mean(x); xm # moyenne des x

[1] 5.5

> ym=mean(c$V2); ym # ordonnée en xm d’apres le cours

[1] 5.514987

> ecart=1.03/sqrt(50) ; ecart*2.01 # demi-largeur de l’IC à 95\% en xm

2
[1] 0.2927846

> ym-2.01*ecart; ym+2.01*ecart # bornes de l’IC à 95\% en xm

[1] 5.222202

[1] 5.807772

13. On se sert des valeurs données par R dans le summary, à la ligne intercept, pour faire le même calcul
que précédemment : β0 ∈ 0.4828 ± 2.01 × 0.3146 = 0.48 ± 0.63. On remarque que l’IC est plus de
deux fois plus large au point x = 0 que x̄, à cause de l’incertitude sur la pente.
14. On remarque que 0 est situé dans l’IC de β0 , donc la valeur numérique obtenue n’est pas significati-
vement différente de 0 avec un risque α = 5%. Cela correspond bien au résultat de summary, où la
p-value associée à H0 : β0 = 0 a une valeur de 0.131. On peut donc envisager un modèle alternatif
y = βx+, qui correspond à cette hypothèse nulle. Le choix entre ces deux modèles a une interprétation
expérimentale. Le modèle passant par 0 correspond à la loi de conductivité la plus simple, où la pente
peut être interprétée comme la conductivité molaire ionique de l’espèce chimique testée. Dans l’autre
modèle, l’ordonnée à l’origine pourrait indiquer un “bruit de fond” expérimental, p. ex. dû à un mauvais
étalonnage de nos instruments (erreur systématique), mais pourrait aussi être un effet réel, indiquant p.
ex. que notre solvant n’a pas une conductivité nulle (p. ex. eau non pure). Cette conductivité résiduelle
n’est pas significativement différente de 0 compte tenu de notre précision expérimentale, mais pas loin,
donc il faudrait éventuellement en tenir compte pour des mesures précises... le choix entre les modèles
n’est donc pas un problème purement statistique.
15. Sans surprise, puisque le paramètres n’était pas significatif, les graphes ne sont pas très différents, et la
légère pente visible sur les résidus n’est sans doute pas elle-même significative.
> c<-read.table("conduct.txt")
> l0=lm(c$V2~c$V1+0); par(mfrow=c(1,2))
> plot(c$V1, c$V2); abline(b=l0$coefficients,a=0)
> plot(c$V1,l0$residuals)
> abline(h=0,lt="dashed")

● ●
10

● ● ●
2

● ● ● ●
● ●

● ● ●
8

● ● ● ●
● ●
● ● ●

l0$residuals


1

● ● ●
● ● ●
● ●
● ●
6
c$V2

● ●
● ● ● ● ●



0

● ● ●
● ●
● ● ●
● ●
4

● ● ● ● ● ● ● ●
● ● ● ●
● ●
● ● ●
−1

● ● ●
● ● ● ●

2

● ● ●
● ● ●
−2



● ●
0

2 4 6 8 10 2 4 6 8 10

c$V1 c$V1

16. On utilise le résultat du cours sur la prédiction d’une mesure en x̄, qui suit une loi de Student avec
V(ȳ1 ) = σ(1 + 1/n) (attention à prendre la racine pour avoir l’écart-type). Ici on utilise l’estimation de
σ. Pour les points différents de x̄, il faut ajouter le terme dû à la pente, qui est un peu plus compliqué à
calculer mais pas tellement grâce à R :
> ecart=1.03*sqrt(1+1/50); 2.01*ecart # demi-largeur au point xbar

3
[1] 2.090901

> ym-2.01*ecart; ym+2.01*ecart

[1] 3.424086

[1] 7.605888

> t=(8.5-xm)**2/sum((x-xm)**2); ordonnee=0.4828+0.9150*8.5; ecart=1.03*sqrt(1+1/50+t)


> ordonnee-2.01*ecart; ordonnee+2.01*ecart # au point x=8.5

[1] 6.147155

[1] 10.37344

> t=(15-xm)**2/sum((x-xm)**2); ordonnee=0.4828+0.9150*15; ecart=1.03*sqrt(1+1/50+t)


> ordonnee-2.01*ecart; ordonnee+2.01*ecart # au point x=15

[1] 11.90354

[1] 16.51206

2 Danger de la régression
Pour aller plus vie, je trace ici successivement les graphes des points avec la droite de régression, et le
graphe des résidus associé.
1. Pas de problème apparent, on peut faire la régression, sachant juste que le faible nombre de résidus rend
délicate la vérification de la normalité.
2. Absence de linéarité, visible dans la forme de parabole des résidus. On le corrigerait par une régression
non linéaire en ajoutant un terme en x2 dans la fonction.
3. Un point aberrant, qui a un effet de levier important sur la droite de régression. Tous les résidus sont
décalés. Il faudrait l’enlever des données.
4. Non-validation de l’homéoscédasticité : l’écart-type semble augmenter avec x. Il est difficile de corriger
ce problème dans le cadre de la modélisation vue en cours.
5. Point aberrant, cf 3.
6. Point aberrant, cf 3.

> par(mfrow=c(4,3),mar=c(4,5,1,0))
> d=read.table(’danger.txt’)
> x=d$V1
> for (i in 1:3) {
+ y=d[,i+1]; plot(x,y); l=lm(y~x); abline(l); text(8,max(y)-2,i)}
> for (i in 1:3) {
+ y=d[,i+1]; l=lm(y~x); plot(x,l$residuals); abline(h=0,lt=’dashed’)}
> for (i in 4:6) {
+ y=d[,i+1]; plot(x,y); l=lm(y~x); abline(l); text(8,max(y)-2,i)}
> for (i in 4:6) {
+ y=d[,i+1]; l=lm(y~x); plot(x,l$residuals); abline(h=0,lt=’dashed’)}

4
30
● ● ● ●

20
● ●
1 2 ● ● 3

15

25
● ● ● ● ●
● ●

15

20
10
● ●
y

y

● ●

15
10

● ● ●
● ● ●

5
● ● ●
● ● ●
● ● ● ●

10
● ●
● ●
5

● ● ●

0
10 15 20 10 15 20 10 15 20

x x x

● ● ●

4

4

● ●

8
● ● ● ●

2
● ●

l$residuals

l$residuals

l$residuals
2

6


0

4
● ●
0

● ● ●
● ●

−6 −4 −2

2

● ●
−4 −2

● ●
● ● ● ●
● ●
● ●
● ● ● ●

−2
● ●

● ● ● ●

10 15 20 10 15 20 10 15 20

x x x
20

● ● ●
● ●
5 6

25
4
30




15

20
25

● ● ●
y

y
10

15
20

● ●

● ●
● ● ● ● ● ●
15

10
5

● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ● ● ● ●
● ● ●

10 15 20 10 15 20 10 15 20

x x x
15

● ● ●

4


10


10


2

● ●
l$residuals

l$residuals

l$residuals
0

● ● ●
5
5
−2

● ● ● ●
● ● ● ●
● ● ● ● ●
● ● ● ●
0
0

● ● ● ●

● ●
−6

● ●
● ● ●
● ● ●
● ●

10 15 20 10 15 20 10 15 20

x x x

3 Test de validation des hypothèses de régression


3.1 Linéarité
3.2 Homéoscédasticité
Hypothèse nulle : la variance (de population) est identique dans les différents groupes (il calcule alors la
probabilité d’observer des variances d’échantillons aussi différentes que celles observées). Le test donne une
p-value de 0.67, donc nous n’avons aucune raison de rejeter l’hypothèse nulle, et le test confirme donc notre
impression visuelle. Remarquez que, puisque les résidus et les points de données diffèrent juste par la valeur de
la droite, les variances des différents groupes sont les mêmes dans les deux cas, et on aurait donc pu appliquer
le test directement sur les points de données.
> c<-read.table("conduct.txt"); l=lm(c$V2~c$V1); bartlett.test(l$residuals~c$V1)

Bartlett test of homogeneity of variances

data: l$residuals by c$V1


Bartlett’s K-squared = 6.6668, df = 9, p-value = 0.6718

3.3 Normalité
L’hypothèse nulle est que l’échantillon testé suit une loi normale. On l’applique sur l’ensemble des résidus.
Contrairement au cas précédent, il ne faut surtout pas l’appliquer sur les points de données, puisque ce sont les

5
résidus qui suivent une loi normale centrée.

> shapiro.test(l$residuals)

Shapiro-Wilk normality test

data: l$residuals
W = 0.96908, p-value = 0.2123

3.4 Test de normalité 2


> par(mfrow=c(1,3),mar=c(4,5,3,0))
> a=runif(10,5,10); qqnorm(a); qqline(a); shapiro.test(a)

Shapiro-Wilk normality test

data: a
W = 0.87394, p-value = 0.1111

> a=runif(30,5,10); qqnorm(a); qqline(a); shapiro.test(a)

Shapiro-Wilk normality test

data: a
W = 0.92454, p-value = 0.03523

> a=runif(100,5,10); qqnorm(a); qqline(a); shapiro.test(a)

Shapiro-Wilk normality test

data: a
W = 0.94707, p-value = 0.0005345

Normal Q−Q Plot Normal Q−Q Plot Normal Q−Q Plot


10

10

● ● ● ●
● ●●●●●●●
● ● ●● ●


●●
9


●●



●●

●●
9

●●●● ●
9


● ● ●
●● ●


Sample Quantiles

Sample Quantiles

Sample Quantiles

● ● ●●


8

● ●


●●● ●




8




●●
8

● ●

● ● ●






●●

7


●●

7

● ●●

7

●●

● ●



●● ●




●●
●●
● ●
6


6

●●

● ●●
●●


6

● ●

● ●●
● ●●
●●●
● ● ● ●
5

−1.5 −0.5 0.5 1.0 1.5 −2 −1 0 1 2 −2 −1 0 1 2

Theoretical Quantiles Theoretical Quantiles Theoretical Quantiles

Voir code pour la réponse aux questions. La p-value obtenue dépend bien sûr de l’échantillon aléatoire tiré,
surtout pour les faibles nombres de tirages. Mais on constate qu’avec dix points, il est absolument impossible de
mettre en évidence la non-normalité. Avec 30 points, c’est limite et dépend beaucoup de l’échantillon tiré (que
ce soit graphiquement ou avec Shapiro). Avec 100 points, on y arrive bien. Ça montre que le test de normalité
nécessite beaucoup de points pour avoir une puissance suffisante, et dans le cas d’une distribution uniforme,
plusieurs dizaines. Pour une distribution qui serait bien plus éloignée encore d’une loi normale (p. ex. une
distribution bimodale), moins de points pourraient suffire.

Vous aimerez peut-être aussi