Sam Meyer - 4BB-BIOSTAT2 2020-2021

TP1 : Régression linéaire 2

1 Régression multiple appliquée au modèle polynomial

> c<-[Link]("[Link]",header=T)
> par(mfrow=c(1,2)); plot(c$x, c$y); l=lm(c$y~c$x)
> abline(l$coefficients); plot(c$x,l$residuals); abline(h=0,lt="dashed")
60

●
●
● ● ●
● ●

6
●
● ●
50

●
● ● ● ●
●
●

4
●
● ●
l$residuals
●
40

● ●
●
2
● ●
●
c$y

●
● ● ●
●
30

● ●
0

● ● ●
●
● ● ●
● ● ● ●
−2

● ● ●
20

●
● ●
●
● ● ●
● ●
● ● ●
● ●
● ● ● ●
●
●
10

● ●
● ●
● ●
●
● ●
−6

2 4 6 8 10 2 4 6 8 10

c$x c$x

1. La répartition des points ne semble pas linéaire.

2. On le confirme avec le graphe des résidus, qui sont répartis suivant une parabole et non autour de 0.
3. Une parabole indique une dépendance en x2 , donc l’ajout d’un terme quadratique pourrait supprimer
cette tendance dans les résidus et corriger le problème.
4. Voir code. Le graphe des résidus (ou ceux fournis par R) indique que les résidus sont bien centrés et de
variance constante. Le graphe des quantiles normaux indique que la normalité ne pose pas de problème.
Donc ce modèle est valable.
5. Le paramètre d’ordre 1 n’est pas significatif. Sauf raison majeure, on préfère un modèle plus simple
avec moins de paramètres. Ici la p-value est de 0.279 pour l’hypothèse nulle où ce paramètre est nul,
et donc avec un risque α de 5%, on privilégie le modèle à deux paramètres sans terme linéaire, et on
recalcule les paramètres.

> c<-[Link]("[Link]",header=T); x=c$x; x2=x*x; par(mfrow=c(1,2)); plot(c$x, c$y);

> l=lm(c$y~x+x2); summary(l); coef=l$coefficients;

Call:
lm(formula = c$y ~ x + x2)

Residuals:
Min 1Q Median 3Q Max
-1.92555 -0.66764 -0.04034 0.60382 2.68722

Coefficients:
Estimate Std. Error t value Pr(>|t|)

1
(Intercept) 7.26603 0.53535 13.572 <2e-16 ***
x 0.24473 0.22359 1.095 0.279
x2 0.47820 0.01981 24.141 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.018 on 47 degrees of freedom

Multiple R-squared: 0.9963, Adjusted R-squared: 0.9961
F-statistic: 6325 on 2 and 47 DF, p-value: < 2.2e-16
> xs=seq(1,10,.1) # on trace la courbe en spécifiant une suite de points x
> lines(xs,l$coef[1]+coef[2]*xs+coef[3]*xs*xs)
> plot(x,l$residuals); abline(h=0,lt="dashed")
> l2=lm(c$y~x2); summary(l2);
Call:
lm(formula = c$y ~ x2)

Residuals:
Min 1Q Median 3Q Max
-1.79609 -0.60305 -0.03271 0.62234 2.55967

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.798512 0.223925 34.83 <2e-16 ***
x2 0.499332 0.004449 112.23 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.02 on 48 degrees of freedom

Multiple R-squared: 0.9962, Adjusted R-squared: 0.9961
F-statistic: 1.26e+04 on 1 and 48 DF, p-value: < 2.2e-16
60

●
●
● ●
●
50

●
●
●
● ●
● ●
● ● ●
l$residuals

● ●
40

● ●
1

● ● ●
● ●
c$y

●
● ● ●
● ● ● ● ●
●
30

● ● ● ●
0

● ● ● ● ● ●
● ●
● ●
● ● ●
● ●
20

●
● ● ● ●
●
−1

● ● ● ●
● ●
● ● ●
● ●
●
10

●
● ●
● ●
●
● ●
−2

2 4 6 8 10 2 4 6 8 10

c$x x

2 Transformation de variable
> par(mfrow=c(1,2),mar=c(4,4,1,0))
> b<-[Link]("[Link]",header=T); t=b$t; n=b$n; nl=log(n)

2
> plot(t,n); plot(t,nl,ylab=’ln n’); [Link](n~t)
Bartlett test of homogeneity of variances

data: n by t
Bartlett’s K-squared = 280.94, df = 9, p-value < 2.2e-16
> l=lm(nl~t); summary(l); coef=l$coefficients; abline(coef)
Call:
lm(formula = nl ~ t)

Residuals:
Min 1Q Median 3Q Max
-0.67627 -0.17986 -0.00648 0.21567 0.57763

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.93074 0.09648 9.647 8.11e-13 ***
t 1.01278 0.01555 65.132 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3158 on 48 degrees of freedom

Multiple R-squared: 0.9888, Adjusted R-squared: 0.9886
F-statistic: 4242 on 1 and 48 DF, p-value: < 2.2e-16
> [Link](nl~t); [Link](l$residuals)
Bartlett test of homogeneity of variances

data: nl by t
Bartlett’s K-squared = 12.758, df = 9, p-value = 0.1739
Shapiro-Wilk normality test

data: l$residuals
W = 0.97901, p-value = 0.5109

● ●
● ●
●
● ●
●
10

●
●
● ● ●
●
●
60000

●
●
8

● ● ●
●
●
ln n

●
n

●
6

●
●
● ●
●
● ●
● ●
20000

●
4

● ●
●
●
● ● ●
●
●
●
●
● ●
2

● ● ●
●
● ● ● ● ● ●
0

2 4 6 8 10 2 4 6 8 10

t t

3
 1. n(t) = n0 exp(ln2t/T ) où n0 est l’inoculum.
2. A première vue, ça pourrait correspondre à une exponentielle.
3. La variance des points augmente avec t, ce qui ne peut être corrigé par la méthode précédente, même
en choisissant une fonction f (t) qui collerait le mieux possible à nos points.
4. Le test de Bartlett sur les données brutes conduit en effet à rejeter complètement l’homéoscédasticité.
5. On travaille sur le logarithme du nombre de bactéries, ce qui devrait régulariser ces données. Cela est
confirmé par les graphes de validation du modèle créés par R, ou par les tests réalisés sur les résidus
(voir code).
6. On calcule d’abord l’IC sur la pente : β1 ∈ 1.013 ± 2.01 ∗ 0.0155 = [0.98, 1.044]. Ici 2.01 est le
résultat de qt(.975,48). On utilise ensuite la question 1 pour passer de la pente à T : avec ln n(t) =
β0 + β1 t + , on a donc n(t) = eβ0 eβ1 t e , et donc T = ln2/β1. Cette fonction est décroissante,
donc si 95% des valeurs de β1 sont comprises dans [0.98,1.044], on peut juste calculer les deux bornes
correspondantes pour T (attention, ça s’inverse) : T ∈ [0.66, 0.71] (en heures). Remarquez que, comme
la fonction inverse n’est pas linéaire, la distribution de T n’est plus symétrique autour de la valeur
centrale.
7. On doit appliquer la formule du cours pour l’intervalle de confiance (pas d’erreur de prédiction) au
point t = 15, comme au TP1. Attention, après application de l’exponentielle, l’IC n’est plus du tout
symétrique ! Cf code :
> # droite passe par le point
> m=[Link](coef[1]+coef[2]*15); m
[1] 16.12247
> # erreur standard? Formule de la variance
> es=sqrt(.3158**2*(1/50+(15-5.5)**2/sum((t-5.5)**2)))
> # IC sur ln n
> m-2.01*es; m+2.01*es
[1] 15.81229
[1] 16.43265
> # IC sur n
> exp(m-2.01*es); exp(m+2.01*es)
[1] 7365310
[1] 13696521

3 Le criquet thermomètre
1. > a<-[Link]("[Link]",header=T)
> plot(a)
> # pas de violation manifeste, mais la distribution des points n’est pas
> # homogène, donc il faut vérifier en plottant les résidus
34

●
30

●
● ●
Temp

●
● ●
●
26

●
●
22

● ●
● ●

15 16 17 18 19 20

Freq

4
2. On voit que les deux variables semblent corrélées, il pourrait donc y avoir une relation linéaire entre les
deux mais il faut le tester. Attention, ici les deux variables sont mesurées dans différentes conditions,
aucune n’est contrôlée.
3.
4. Sachant qu’on va vouloir déterminer la température en fonction de la fréquence mesurée, on considère
cette dernière comme variable explicative. Mais attention, la condition de “variable contrôlée” n’est pas
vérifiée ! Dans la suite, on va supposer que si la mesure de la fréquence est précise, le fait de ne pas
contrôler la fréquence n’empêche pas l’application de la régression. Mais attention, on devra toujours
garder ce problème en tête ! L’énoncé précise qu’on a pris des précautions pour assurer des mesures
indépendantes, là-dessus puisqu’on n’a pas d’autres informations, on est obligé de faire confiance à
l’expérimentateur.
Pour les autres hypothèses, on va regarder les graphes (i) des résidus (résidus centrés en 0), (ii) quan-
tile normal-quantile (normalité), (iii) des résidus ou le troisième graphe sur l’amplitude des résidus
(homéoscédasticité). Aucune des hypothèses ne semble violée, compte tenu du faible nombre de points.
> a<-[Link]("[Link]",header=T)
> x=a$Freq
> y=a$Temp
> par(mfrow = c(2, 2))
> N=length(x)
> lm1=lm(y~x)
> plot(lm1)
> # Propriétés des résidus
> summary(lm1)
Call:
lm(formula = y ~ x)

Residuals:
Min 1Q Median 3Q Max
-3.3123 -1.5979 -0.2819 1.6998 3.1057

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.0928 5.8319 -0.873 0.398364
x 1.8951 0.3503 5.410 0.000119 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.245 on 13 degrees of freedom

Multiple R-squared: 0.6925, Adjusted R-squared: 0.6688
F-statistic: 29.27 on 1 and 13 DF, p-value: 0.000119

5
 Residuals vs Fitted Normal Q−Q

1.5
● 11 11 ●
●

Standardized residuals
● ● ●
2
●
● ● ● ●
●

0.5
Residuals

●
●
0

● ●
● ● ●

−0.5
●
●
●
−2

●
●
●
●
●9 2●

−1.5
●9 ●2
−4

22 24 26 28 30 32 −1 0 1

Fitted values Theoretical Quantiles

Scale−Location Residuals vs Leverage

1
●9 2●
1.2

● 11

1.5
● 11
● ●
Standardized residuals

● ● 0.5
Standardized residuals
●

● ● ● 3●
● ●
●
0.5
0.8

●
●

● ● ●
●
−0.5

●
0.4

●
●
● 1●

●
−1.5

0.5
● ●
Cook's distance
0.0

22 24 26 28 30 32 0.00 0.10 0.20 0.30

Fitted values Leverage

5. On teste l’hypothèse nulle β = 0, la p-value donnée par R suite au test de Student associé (portant
sur la variable aléatoire b) est d’environ 10−4 donc on peut rejeter l’hypothèse nulle avec un risque de
première espèce faible. La température a un effet significiatif sur la fréquence.
6. La précision est maximale à la moyenne des xi . D’après le cours,√la température moyenne à cette fré-
quence vaut ȳ et l’écart-type de cette variable aléatoire est de σ/ N . L’intervalle de prédiction n’est
pas le même, puisqu’il inclut notre erreur de calibration ET le bruit inhérent à une mesure individealle.
On utilise la formule vue en cours, ou la commande de R qui donne directement le même résultat.
Comme vous le notez, l’erreur est considérablement plus grande. Remarquez que j’ai fait l’approxima-
tion que l’intervalle à 95% est deux fois l’écart-type de chaque côté. Plus précisément, on peut calculer
sa largeur en étudiant la distribution de Student à N − 2 ddl : qt(0.975,df=13)=2.16 au lieu de
2 écarts-type. La définition de ces deux intervalles : intervalles où on estime à 95% la probabilité de
trouver le paramètre α et la “vraie” température lors d’une mesure individuelle, respectivement.
> # 2. La précision est maximale à la moyenne des xi:
> xbar<-mean(x)
> # ecart type des erreurs
> sigma=2.245
> # Température moyenne correspondant à 14 Hz
> # et intervalle de confiance
> # largeur de part et d’autre
> mean(y)-2*sigma/sqrt(N); mean(y)+2*sigma/sqrt(N)
[1] 25.14402
[1] 27.46265
> # Intervalle de prédiction
> mean(y)-2*sigma*sqrt(1+1/N) ; mean(y)+2*sigma*sqrt(1+1/N)
[1] 21.66608
[1] 30.94059

6
 7. A 14 Hz, on reprend la même formule
> c=lm1$coef # coefficients du modèle $
> me=[Link](c[1]+14*c[2]) # valeur moyenne de la température à 14Hz
> me
[1] 21.43914
> # pour l’incertitude, on a besoin de la somme des (xi-xbar)**2
> sx=sum((x-xbar)**2)
> # on calcule le double de l’écart-type, par la formule de la variance
> er=2*sigma*sqrt(1+1/N+(14-xbar)**2/sx)
> # finalement l’intervalle est donné par
> me-er;me+er
[1] 16.46545
[1] 26.41283
> # Pour la question 8, on veut calculer la précision avec 5 criquets: on utilise
> # donc la formule du cours avec une valeur m=5:
> er=2*sigma*sqrt(1/5+1/N+(14-xbar)**2/sx)
> me-er;me+er
[1] 18.50495
[1] 24.37334
> # déjà bien plus précise !
8. On fait ce choix (i) pour gagner en précision car on veut faire la moyenne de 5 mesures plutôt que
d’utiliser une mesure individuelle ; (ii) pour assurer l’indépendance des résultats, plutot que de faires
des mesures successives sur un même criquet. Remarquez que cette indépendance est toute relative, si
en réalité la réponse des criquets dépend de l’heure de la journée par exemple ! ! ! Cf résultat ci-dessus,
la précision est effectivement bien meilleure (facteur plus que 2).
9. Il faut bien réfléchir à la température que nous cherchons : ce n’est pas la température où le criquet
chante en moyenne à 10 Hz ! ! ! Sinon le risque de bruitage est de 1/2, car il a une chance sur deux
de chanter à une fréquence audible (si il est un peu au-dessus plutôt qu’en dessous). Il faut trouver la
température maximale où le criquet a 95% de chances de chanter en-dessous de 10 Hz. On va donc
tracer, pour chaque fréquence, l’intervalle de prédiction à 90% : par définition, cela correspond à deux
courbes telles qu’on ait 5% des valeurs au-dessus, et 5% des valeurs en-dessous. Ces courbes sont sur
le graphe suivant, j’ai utilisé la formule du cours. Remarquez qu’il ne s’agit pas de droites, mais de
branches de parabole, on le voit beaucoup mieux pour l’intervalle sur les paramètres du modèle car leur
estimation est beaucoup plus précise dans la zone de calibration.
Graphiquement, la température recherchée est celle qui coupe la branche de parabole inférieure à une
fréquence de 10 Hz. En effet, on s’attend alors à ce que la partie gauche (non audible) contienne bien
95% des valeurs de fréquence à cette température, on aura donc absence de bruitage avec un risque de
5%. On trouve une valeur de 8.1°C. Toutefois, ce calcul n’est pas rigoureux puisque les intervalles de
confiance ne portent que sur la variable dépendante : ici on cherche en fait un intervalle de confiance sur
la fréquence, donc en toute rigueur il faudrait refaire la régression en prenant cette fois la température
comme variable explicative, et chercher pour quelle température on a la branche inférieure qui passe à
la fréquence de 10 Hz... Je vous en fais grâce.
10. Plus difficile : pour que les 5 criquets soient silencieux avec 95% de confiance (risque α = 5%), chacun
individuellement doit être silencieux avec environ 99% de confiance ! (cf Bonferroni). On va donc cette
fois considérer la courbe rose, telle que 1% des valeurs soient en-dessous (IC à 98% sur les prédictions.
Elle passe à 10 Hz pour une température de 5°C (ligne rose). Le thermomètre à 5 criquets n’est donc
vendable que jusqu’à cette température, mais il a l’avantage d’être plus précis que celui à un criquet (cf
ci-dessus) ! Il est vrai que (i) un thermomètre ne marchant que jusqu’à 5°C n’a pas beaucoup d’intérêt

7
pour l’homme et (ii) sans être un spécialiste des insectes, je doute de leur longévité à cette température.
En conclusion, tout cela n’a pas un grand intérêt pratique, comme vous pouviez vous en douter.
Petites remarques importantes sur tous ces graphes : (1) Au centre du graphe, l’intervalle de confiance
sur la réponse moyenne est bien plus réduit que la dispersion des points initiaux : c’est tout l’intérêt de
combiner plusieurs mesures ! (2) Loin des points, le modèle devient vite imprécis. (3) Pour la prédiction
du résultat d’une seule mesure, au centre c’est la variabilité σ qui domine (le rouge est beaucoup plus
large que le bleu), en revanche sur les côtés l’imprécision du modèle est au moins aussi grande (les
courbes sont plus proches).
> a<-[Link]("[Link]",header=T)
> x=a$Freq
> y=a$Temp
> sigma=2.245
> # v=sigma**2
> xbar=mean(x)
> lm1=lm(y~x)
> c=lm1$coef
> sx=sum((x-xbar)**2)
> # $on va tracer les intervalles pour des valeurs de x comprises entre 9 et 22 Hz
> # je crée un tableau avec ces valeurs de x, séparées de 0.5
> plotx <- seq(5,max(x)+2,.5)
> n=length(plotx)
> # je calcule les intervalles de prédiction correspondants et la temp max
> pred1=[Link](c[1]+c[2]*plotx)-qt(.95,df=13)*sigma*sqrt(1+1/15+(plotx-xbar)**2/sx)
> pred2=[Link](c[1]+c[2]*plotx)+qt(.95,13)*sigma*sqrt(1+1/15+(plotx-xbar)**2/sx)
> tmax=[Link](c[1]+c[2]*10)-qt(.95,13)*sigma*sqrt(1+1/15+(10-xbar)**2/sx)
> tmax
[1] 8.074565
> # pour la boite à 5 criquets, il faut un intervalle à 99\%
> # compte tenu de la multiplicité des criquets
> pred51=[Link](c[1]+c[2]*plotx)-qt(.99,df=13)*sigma*sqrt(1+1/15+(plotx-xbar)**2/sx)
> pred52=[Link](c[1]+c[2]*plotx)+qt(.99,13)*sigma*sqrt(1+1/15+(plotx-xbar)**2/sx)
> # pour le plaisir, voici l’intervalle de confiance sur le modèle
> conf1=[Link](c[1]+c[2]*plotx)-qt(.95,df=13)*sigma*sqrt(1/15+(plotx-xbar)**2/sx)
> conf2=[Link](c[1]+c[2]*plotx)+qt(.95,13)*sigma*sqrt(1/15+(plotx-xbar)**2/sx)
> # et je trace tout
> plot(x,y,style=’p’,xlim=c(8,22),ylim=c(5,35))
> lines(plotx,pred1,col=4)
> lines(plotx,pred2,col=4)
> lines(plotx,conf1,col=3)
> lines(plotx,conf2,col=3)
> lines(plotx,pred51,col=6)
> lines(plotx,pred52,col=6)
> abline(v=10,col=1)
> abline(h=tmax,col=4,lty=1)
> abline(h=5,col=6,lty=1)

8
 35
●

●
30

●
● ●
●
● ●
●
25

●
●

●●
● ●
20
y

15
10
5

8 10 12 14 16 18 20 22