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

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

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