Statistique théorique et appliquée - Tome 2

 
Utilisation de R pour les exemples 14.3.1 à 14.3.3

Exemple 14.3.1
Exemple 14.3.2
Exemple 14.3.3

par Emmanuel Nowak

 

Exemple 14.3.1. Décroissance des restes de fongicide observés sur laitue : ajustement d'une droite de régression

Il s'agit ici d'étudier les restes de fongicide observés sur laitue en fonction du nombre de jours après application. Il faut tout d'abord importer les données, qui sont celles relatives à l'exemple 3.6.1 (fichier s2e03061), et diviser la valeur des restes par 100 afin d'obtenir les données réelles :

s2e03061 <- read.table("C:/Dagnelie/st2donn/txt.2/s2e03061.txt",sep="\t",header=T)
s2e03061$Reste <- s2e03061$Reste/100
s2e03061[1:5,]

  Jours Reste
1     1  8.96
2     1 11.50
3     1 13.12
4     3  6.72
5     3  7.68 

On peut alors construire le modèle de régression sur les données transformées (transformation logarithmique pour les restes) :

attach(s2e03061)
Z <- log10(Reste)
model.Z <- lm(Z~Jours)
resum <- summary(model.Z)
resum

Call:
lm(formula = Z ~ Jours)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.105106 -0.054282 -0.003152  0.065187  0.116586 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.053333   0.037356   28.20 7.32e-11 ***
Jours       -0.051985   0.004679  -11.11 6.00e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.08053 on 10 degrees of freedom
Multiple R-Squared: 0.9251,     Adjusted R-squared: 0.9176 
F-statistic: 123.5 on 1 and 10 DF,  p-value: 6.004e-07

La représentation graphique des données (figure 3.6.1 du livre), avec ou sans transformation, ainsi que la droite de régression, s'obtiennent de la fa┴on suivante :

par(mfrow=c(1,2))
plot(Reste~Jours,ylab="Fongicide (p.p.m.)")
plot(Z~Jours,ylab="Log10(fongicide)")
abline(model.Z)

s2e03061.jpg (45,5 Ko)

Toutes les informations nécessaires à la poursuite de l'exemple figurent dans la sortie précédente. Ces valeurs sont contenues dans l'objet 'resum' dont on peut lister le contenu :

names(resum)

[1] "call"          "terms"         "residuals"     "coefficients" 
[5] "aliased"       "sigma"         "df"            "r.squared"    
[9] "adj.r.squared" "fstatistic"    "cov.unscaled" 

Tout d'abord, 'sigma' contient l'estimation de l'écart-type résiduel. La variance résiduelle est donc estimée par  :

resum$sigma^2

[1] 0.006484994

Afin d'en calculer les limites de confiance, il faut déterminer les valeurs relatives à une variable khi-carré dont le nombre de degré de liberté est celui de la variance résiduelle, c'est-à-dire 10 comme indiqué plus haut, pour un degré de confiance égal à 95% :

chi2.alpha <- c(qchisq(0.025,df=10,lower.tail=F),qchisq(0.025,df=10,lower.tail=T))
chi2.alpha

[1] 20.483177  3.246973

Il ne reste plus qu'à multiplier l'estimation de la variance résiduelle par son nombre de degré de liberté et diviser par les valeurs obtenues ci-dessus pour obtenir les limites de confiance :

10*resum$sigma^2/chi2.alpha

[1] 0.00316601 0.01997243

Les estimations de l'ordonnée à l'origine et de la pente de la droite de régression ainsi que leurs écart-types peuvent être récupérés afin de déterminer les limites de confiance de ces coefficients :

resum$coef

               Estimate  Std. Error   t value     Pr(>|t|)
(Intercept)  1.05333259 0.037356459  28.19680 7.315106e-11
Jours       -0.05198501 0.004678704 -11.11098 6.004335e-07

Il faut alors déterminer les valeurs relatives à une variable t de Student à 10 degrés de liberté, pour un degré de confiance égal à 95%, et appliquer la formule habituelle :

t.alpha <- qt(0.025,df=10,lower.tail=F)
resum$coef[,1]-t.alpha*resum$coef[,2]; resum$coef[,1]+t.alpha*resum$coef[,2]

(Intercept)       Jours 
 0.97009721 -0.06240981 
(Intercept)       Jours 
 1.13656797 -0.04156021
 
Haut de la page




Exemple 14.3.2. Décroissance des restes de fongicide observés sur laitue : étude des résidus

Le logiciel R dispose de nombreux outils diagnostiques mais la terminologie utilisée ne correspond pas toujours à celle du livre. En reprenant les données de l'exemple 14.3.1, deux modèles peuvent être construits, l'un sur les données initiales et l'autre sur les données transformées :

model.Y <- lm(Reste~Jours)
model.Z <- lm(Z~Jours)

Les valeurs hi, arrondies ci-dessous à quatre décimales, sont les éléments diagonaux de la 'hat' matrice et sont donc les mêmes pour les deux modèles :

round(hatvalues(model.Y),digit=4) # ou hatvalues(model.Z)

     1      2      3      4      5      6      7      8      9     10     11     12 
0.1764 0.1764 0.1764 0.1190 0.1190 0.1190 0.0852 0.0852 0.0852 0.2861 0.2861 0.2861 

Les résidus di, arrondis à deux ou trois décimales, peuvent être obtenus directement pour chacun des deux modèles :

round(model.Y$res,digit=2)

    1     2     3     4     5     6     7     8     9    10    11    12 
-0.72  1.82  3.44 -1.70 -0.74 -0.42 -2.06 -2.06  0.18  0.43  0.75  1.07
round(model.Z$res,digit=3)

     1      2      3      4      5      6      7      8      9     10     11     12 
-0.049  0.059  0.117 -0.070 -0.012  0.006 -0.105 -0.105  0.094 -0.042  0.025  0.083 

Les résidus di", nommés résidus studentisés dans le livre, sont appelés résidus standardisés dans R :

round(rstandard(model.Y),digit=2)

    1     2     3     4     5     6     7     8     9    10    11    12 
-0.46  1.16  2.20 -1.05 -0.46 -0.26 -1.25 -1.25  0.11  0.30  0.52  0.74 
round(rstandard(model.Z),digit=2)

    1     2     3     4     5     6     7     8     9    10    11    12 
-0.67  0.81  1.60 -0.93 -0.16  0.08 -1.36 -1.36  1.23 -0.62  0.36  1.22

Les représentations graphiques de ces résidus en fonction du temps sont les suivantes :

plot(rstandard(model.Y)~Jours,ylab="Résidus d'' de y")
abline(h=0,lty=2)
plot(rstandard(model.Z)~Jours,ylab="Résidus d'' de z")
abline(h=0,lty=2)

s2e14031.jpg (42,4 Ko)

Notons que la commande 'rstudent' de R fournit les résidus studentisés jacknife (voir paragraphe 14.3.3.7° du livre), et que d'autres outils sont également disponibles : consulter la fonction 'influence.measures' pour plus de détails.



Exemple 14.3.3. Décroissance des restes de fongicide observés sur laitue : distances de Cook

Les distances de Cook relatives aux deux modèles précédents, arrondies à deux décimales, s'obtiennent immédiatement :

round(cooks.distance(model.Y),digit=2)

   1    2    3    4    5    6    7    8    9   10   11   12 
0.02 0.15 0.52 0.07 0.01 0.00 0.07 0.07 0.00 0.02 0.05 0.11
round(cooks.distance(model.Z),digit=2)

   1    2    3    4    5    6    7    8    9   10   11   12 
0.05 0.07 0.27 0.06 0.00 0.00 0.09 0.09 0.07 0.08 0.03 0.30

Avant de traiter d'autres exemples, ne pas oublier de détacher le fichier de travail :

detach(s2e03061)

Haut de la page
Autres illustrations avec R




Dernière mise à jour : septembre 2007