Universit Rennes 2 e Statistique M1-GEO

Corrig des TD e
Ouvrages recommands e
Tous ces livres sont ` la BU. Pour les acheter, venir au bureau A-240 ou envoyer un mail : [Link]@[Link] a 1. Agn`s Hamon et Nicolas Jegou, Statistique descriptive. Cours et exercices corrigs., PUR, 2008 e e Pour revoir la base et sinitier ` Rcmdr. a 2. Jrme Pag`s, Statistiques gnrales pour utilisateurs. 1-Mthodologie, PUR, 2005 eo e e e e Transcription du cours donn ` Agrocampus Rennes. Estimation, analyse de variance et rgression puis ea e introduction aux plans dexprience et ` lACP. Introduction ` la statistique pratique, tr`s pdagogique e a a e e et tr`s bien crit. e e 3. Franois Husson et Jrme Pag`s, Statistiques gnrales pour utilisateurs. 2-Exercices et corrigs, c eo e e e e PUR, 2005 Exercices et corrigs en lien avec louvrage prcdent. Quelques TP sur R proposs. e e e e 4. [Link] et al., Statistiques avec R., PUR, 2008 Prsentation du logiciel : objets, graphiques, programmation. Quinze mthodes statistiques classiques e e prsentes avec R. Indispensable pour laspect logiciel. e e 5. [Link]`s, [Link], S.L, Analyse de donnes avec R, PUR 2009. Utile pour la seconde partie du S2 e e e et le master 2.

Rappels

Exercice 1 : Modelisation
Dans cet exercice, nous ralisons certaines lignes de commandes qui taient obtenues directement avec Rcmdr e e au S1. Les commandes de lexercice sont les suivantes : > > > > > > > > > > > #1) a <- 1 b <- 2 n <- 100 X <- seq(0,1,length=n) eps <- rnorm(n,mean=0,sd=0.3) Y <-a*X+b+eps plot(X,Y) model <- lm(Y~X) predict(model) lines(X,predict(model))

On obtient la gure 1. Dans cet exercice, les donnes sont simules selon un mod`le connu qui est ici linaire. En pratique, ce nest e e e e pas le cas : on dispose de donnes mais la relation entre X et Y est inconnue. Comme il ny a que deux e variables en prsence, il est possible de reprsenter simplement les donnes. On peut donc voir facilement si e e e ajuster les points par une droite constitue une simplication en accord avec cette reprsentation. e Si tel est le cas, on postule une relation de la forme Y = aX + b + . (1)

Y 1.5 2.0

2.5

3.0

0.0

0.2

0.4 X

0.6

0.8

1.0

Fig. 1 Ajustement dun mod`le linaire sur des donnes simules. e e e e Deux remarques importantes : En pratique, on ne sait pas si la relation entre X et Y est bien de la forme (1) mais cette formulation prsente le double avantage dtre dinterprtation simple et de conduire ` des calculs ralisables (cf. cours e e e a e du S1). De plus, les param`tres a et b du mod`le (1) sont en pratique bien sr inconnus mais on utilise les donnes e e u e pour les estimer. Nous obtenons une estimation a et des param`tres du mod`le ainsi : b e e > summary(model) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.01196 0.06854 29.355 < 2e-16 *** X 0.93228 0.11842 7.873 4.73e-12 ***

Exercice 2 : Calcul du R2
Pour obtenir la gure 2 nous crivons : e > plot(X,Y) > lines(X,predict(model)) > abline(h=mean(Y),col="red") La valeur de R2 donne par le logiciel est la suivante : e > summary(model) Multiple R-squared: 0.3874, Adjusted R-squared: 0.3812 F-statistic: 61.98 on 1 and 98 DF, p-value: 4.731e-12

Y 1.5 2.0

2.5

3.0

0.0

0.2

0.4 X

0.6

0.8

1.0

Fig. 2 Aide ` linterprtation du R2 . a e Nous la retrouvons en reprenant la formule R2 =

Pn y 2 i=1 y Pn (i )2 (yi ) y i=1

propose : e

> Rdeux <- sum((predict(model)-mean(Y))^2)/sum((Y-mean(Y))^2) > Rdeux [1] 0.3874345 Au terme 1/n pr`s, le dnominateur de R2 correspond ` la variance des yi . Cest donc une mesure de la e e a variabilit totale observe sur les valeurs de Y . Graphiquement, cette quantit sinterpr`te comme la somme e e e e des carrs des longueurs verticales entre les points du nuage et leur projet sur la droite reprsentant la e e e moyenne. Le numrateur lui reprsente la somme des carrs des distances entre les points ajusts par le e e e e mod`le et cette mme droite horizontale : cest une mesure de la variabilit des donnes ajustes. e e e e e On peut montrer que lon a toujours R2 [0, 1] ce qui autorise son interprtation en termes de pourcentages. e Ici, on dira que le mod`le ajust explique 38.7% de la variabilit des donnes. Une valeur proche de 1 signie e e e e que numrateur et dnominateur sont proches : la variabilit des valeurs ajustes est donc proche de celle e e e e des donnes. On conclue alors que le mod`le explique bien la variabilit des donnes. A linverse, une valeur e e e e proche de 0 signie que le dnominateur est bien plus grand que le numrateur. La variabilit des valeurs e e e ajustes est loigne de la variabilit relle des donnes : le mod`le ajust nest pas satisfaisant. e e e e e e e e Remarque 1 Sans insister sur la nature du test eectu, nous voyons que R rend une probabilit critique e e concernant la statistique R2 : p-value: 4.731e-12 Cette probabilit correspond au test de signicativit (ou de nullit) de R2 . Plus prcisment, on teste e e e e e lhypoth`se e H0 : R2 = 0 contre H1 : R2 = 0. Ainsi, on peut considrer que la valeur de R2 observe (0.3874), si elle nest pas vritablement proche de 1, e e e est quand mme tr`s signicativement dirente de 0. e e e

Exercice 3 : Application
Le chier comporte 13 variables. La variable ` expliquer est maxO3 et les autres sont des variables explicatives a potentielles. On dispose de 112 mesures simultanes de ces 13 variables. e 1. On dnit un mod`le linaire entre maxO3 et T12. Nous obtenons la gure 3 et les sorties suivantes : e e e > RegModel.1 <- lm(maxO3~T12, data=Dataset) > summary(RegModel.1) Call: lm(formula = maxO3 ~ T12, data = Dataset) Residuals: Min 1Q -38.0789 -12.7352

Median 0.2567

3Q 11.0029

Max 44.6714

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -27.4196 9.0335 -3.035 0.003 ** T12 5.4687 0.4125 13.258 <2e-16 *** --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 17.57 on 110 degrees of freedom Multiple R-squared: 0.6151,Adjusted R-squared: 0.6116 F-statistic: 175.8 on 1 and 110 DF, p-value: < 2.2e-16 Les param`tres estims sont a = 5.4687 et = 27.4196. e e b 2. De mme avec la variable Ne12 en gure 4 et ci-dessous : e > RegModel.2 <- lm(maxO3~Ne12, data=Dataset) > summary(RegModel.2) Call: lm(formula = maxO3 ~ Ne12, data = Dataset) Residuals: Min 1Q -46.020 -14.487

Median -5.115

3Q 12.571

Max 66.470

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 130.0201 4.9807 26.105 < 2e-16 *** Ne12 -7.9150 0.9042 -8.753 2.77e-14 *** --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 21.74 on 110 degrees of freedom Multiple R-squared: 0.4106,Adjusted R-squared: 0.4052 F-statistic: 76.62 on 1 and 110 DF, p-value: 2.769e-14

maxO3

40

60

80

100

120

140

160

15

20 T12

25

30

Fig. 3 Inuence de T12. Les param`tres estims sont a = 7.9150 et = 130.0201. Le signe de la pente montre que si la e e b nbulosit augmente, la concentration dozone dans lair a tendance ` diminuer. e e a 3. Nous observons que la variable T12 explique mieux la variabilit de maxO3 que Ne12 puisquelle e 2 plus proche de 1 : cest cette variable que lon retiendrait sil fallait choisir. prsente un R e Remarque 2 Nous avons test la signicativit de la pente dun mod`le de rgression simple au e e e e premier semestre. Dans un mod`le de rgression de la forme (1), cela revient a tester e e ` H0 : a = 0 contre H1 : a = 0. Remarquons que la probabilit critique associe a ce test est la mme que celle associe au test de e e ` e e signicativit de R2 . Ainsi, en rgression simple, tester la nullit de R2 quivant a tester la nullit de e e e e ` e a.

Exercice 4 : Analyse de variance

Limportation des donnes montre que la prsence de 3 variables : le poids non-viscr dun poulpe [Link], e e e ee son poids viscr [Link] et son sexe Sexe. La variable ` expliquer est le poids viscr et Sexe est la variable e ee a e ee explicative. Le rsum des donnes montre que cette derni`re est interprte, du fait de son codage, comme e e e e ee une variable quantitative : > summary(Dataset) [Link] [Link] Min. : 40.0 Min. : 20.0 1st Qu.: 350.0 1st Qu.: 300.0 Median : 550.0 Median : 500.0 Mean : 742.3 Mean : 639.5 3rd Qu.: 900.0 3rd Qu.: 800.0 Sexe Min. :1.000 1st Qu.:1.000 Median :2.000 Mean :1.510 3rd Qu.:2.000

maxO3

40 0

60

80

100

120

140

160

4 Ne12

Fig. 4 Inuence de Ne12. Max. :5400.0 Max. :5100.0 Max. :2.000

Il convient de convertir cette variable en facteur : > Dataset$Sexe <- factor(Dataset$Sexe, labels=c(F,M)) > summary(Dataset) [Link] [Link] Min. : 40.0 Min. : 20.0 1st Qu.: 350.0 1st Qu.: 300.0 Median : 550.0 Median : 500.0 Mean : 742.3 Mean : 639.5 3rd Qu.: 900.0 3rd Qu.: 800.0 Max. :5400.0 Max. :5100.0

Sexe F:240 M:250

Nous disposons dun chantillon de 490 poulpes compos de 240 femmelles et 250 mles. Nous pouvons e e a reprsenter la distribution de Y dans chaque modalit de sexe (cf. gure 5) : e e > boxplot([Link]~Sexe, ylab="[Link]", xlab="Sexe", data=Dataset) Nous pouvons calculer les rsums par groupe : e e > numSummary(Dataset[,"[Link]"], groups=Dataset$Sexe, statistics=c("mean", + "sd"), quantiles=c(0,.25,.5,.75,1)) mean sd n F 543.3958 384.0462 240 M 731.8400 722.2932 250

[Link]

1000

2000

3000

4000

5000

F Sexe

Fig. 5 Poids viscr des poulpes selon le sexe. e ee Nous notons M et F le poids moyen respectivement dun poulpe mle et dun poulpe femmelle dans a lensemble des eaux Mauritaniennes. Ces valeurs sont inconnues mais lchantillon nous en donne une e estimation : M = 731.8400 et F = 543.3958. Sur notre chantillon, nous observons donc une dirence qui semble assez sensible entre ces deux valeurs. e e Cette dirence rv`le-t-elle que lon doive considrr que M = F ? Cest la question (entre autres) ` e e e ee a laquelle on peut rpondre en eectuant une analyse de variance avec R. Eectuons cette analyse : e > AnovaModel.1 <- aov([Link] ~ Sexe, data=Dataset) > summary(AnovaModel.1) Df Sum Sq Mean Sq F value Pr(>F) Sexe 1 4348311 4348311 12.848 0.0003717 *** Residuals 488 165155611 338434 --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 > numSummary(Dataset$[Link] , groups=Dataset$Sexe, statistics=c("mean", "sd")) mean sd n F 543.3958 384.0462 240 M 731.8400 722.2932 250 Nous concluons que les dirences de moyennes observes M et F sont signicativement direntes puisque e e e la probabilit critique ` associe au test propos est faible (0.0003717). Au seuil classique de 5%, nous rejetons e a e e H0 : M = F et considrons que le sexe inuence le poids viscr du poulpe. e e ee

Exercice 5 : Approfondissement
En faisant tourner plusieurs fois le programme suivant, on constate que la proportion des carts de moyennes e au moins aussi grands que celui constat ` lexercice prcdent correspond bien ` la probabilit critique e a e e a e trouve : e > sigma.M <- 722.2932 > sigma.F <- 384.0462 > moy.M <- 731.8400 > moy.F <- 543.3958 > [Link] <- abs(moy.M-moy.F) > mu <- 0.5*(moy.M+moy.F) > mu [1] 637.6179 > > > > N <- 100000 > n.M <- 250 > n.F <- 240 > > [Link] <- rep(0,N) > [Link] <- rep(0,N) > [Link] <- rep(0,N) > for(i in 1:N){ + [Link][i] <- mean(rnorm(n.M,mu,sigma.M)) + [Link][i] <- mean(rnorm(n.F,mu,sigma.F)) + [Link][i] <- abs([Link][i]-[Link][i]) + } > > (1/N)*table([Link]>[Link]) FALSE TRUE 0.99962 0.00038

Rgression multivarie e e
1. On crit un mod`le de rgression multiple de la forme : e e e maxO3i = 0 + 1 T 12i + 2 N e12i + i . 2. On peut reprsenter la situation en 3D ainsi : e > scatter3d(Dataset$Ne12, Dataset$maxO3, Dataset$T12, fit="linear", + residuals=TRUE, bg="white", [Link]=TRUE, grid=TRUE, ellipsoid=FALSE, + xlab="Ne12", ylab="maxO3", zlab="T12") La reprsentation est proche de celle faite en gure 6. e Remarque 3 Il est plus dicile en 3D quen 2D de savoir si un mod`le linaire est rellement pertie e e nent. Par ailleurs, d`s que lon a plus de trois variables, toute reprsentation simple devient impossible. e e Le mod`le linaire reste cependant un outil tr`s courant car son interprtation reste assez simple. e e e e

Exercice 6 : Mod`le bivari e e

maxO3

Ne12 T12

Fig. 6 Nuage des points sur trois variables quantitatives. 3. On obtient lquation du plan le plus proche au sens des moindres carrs en ajustant un mod`le linaire e e e e aux observations selon ce crit`re et en lisant les valeurs estimes des param`tres : e e e > LinearModel.1 <- lm(maxO3 ~ T12 +Ne12, data=Dataset) > summary(LinearModel.1) Call: lm(formula = maxO3 ~ T12 + Ne12, data = Dataset) Residuals: Min 1Q -41.13325 -11.79904

Median -0.02899

3Q 9.14797

Max 45.58385

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 7.7077 15.0884 0.511 0.61050 T12 4.4649 0.5321 8.392 1.92e-13 *** Ne12 -2.6940 0.9426 -2.858 0.00511 ** --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 17.02 on 109 degrees of freedom Multiple R-squared: 0.6419,Adjusted R-squared: 0.6353 F-statistic: 97.69 on 2 and 109 DF, p-value: < 2.2e-16 Nous lisons ainsi 0 = 7.7077, 1 = 4.4649 et 2 = 2.6940, param`tres estims qui dnissent e e e

lquation du plan ajust que nous notons : e e Y = 7.7077 + 4.4649T 12 2.6940N e12 4. Les param`tres estims en rgression multiple ne sont pas ceux quon obtiendrait en faisant plusieurs e e e rgressions simples. Il convient dajouter une remarque importante concernant le R2 . e Remarque 4 La valeur de R2 calcule a linstant est suprieure a chacune des deux que nous avions e ` e ` calcul en exercice dapplication en section 1. En fait, il faut savoir que le R2 augmentera toujours e quand on ajoute des variables explicatives. 5. On obtient le graphe des rsidus (cf. gure 7) avec la commande : e > plot(LinearModel.1$residuals)

LinearModel.1$residuals

40 0

20

20

40

20

40

60 Index

80

100

Fig. 7 Reprsentations des rsidus. e e Ce graphe ne rv`le pas de point aberrant ni de structure particuli`re. e e e 6. Pour eectuer une prvisions, nous utilisons les coecients estims : e e > LinearModel.1$coef (Intercept) T12 Ne12 7.707653 4.464851 -2.693975 > LinearModel.1$coef[1]+20*LinearModel.1$coef[2]+4*LinearModel.1$coef[3] (Intercept) 86.22878

Exercice 7 : Avec 3 variables explicatives

1. Le mod`le scrit cette fois : e e maxO3i = 0 + 1 T 12i + 2 N e12i + 3 V x12i + i .

2. On ajuste un mod`le linaire de la mani`re suivante : e e e > LinearModel.2 <- lm(maxO3 ~ T12 +Ne12 +Vx12, data=Dataset) > summary(LinearModel.2) Call: lm(formula = maxO3 ~ T12 + Ne12 + Vx12, data = Dataset) Residuals: Min 1Q -37.4624 -11.4484

Median -0.7219

3Q 8.9082

Max 46.3314

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.8958 14.8243 0.263 0.7932 T12 4.5132 0.5203 8.674 4.71e-14 *** Ne12 -1.6189 1.0181 -1.590 0.1147 Vx12 1.6290 0.6571 2.479 0.0147 * --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 16.63 on 108 degrees of freedom Multiple R-squared: 0.6612,Adjusted R-squared: 0.6518 F-statistic: 70.25 on 3 and 108 DF, p-value: < 2.2e-16 3. On a R2 = 0.6612, valeur suprieure au cas prcdent conformment ` ce qui a t dit en remarque 4. e e e e a ee 4. Lanalyse des rsidus ne rv`le rien de particulier (cf. gure 8) : e e e > plot(LinearModel.2$residuals) 5. Pour eectuer une prvision, nous reprenons les coecients e > LinearModel.2$coef (Intercept) T12 Ne12 Vx12 3.895777 4.513237 -1.618885 1.629015 > LinearModel.2$coef[1]+20*LinearModel.2$coef[2]+ 4*LinearModel.2$coef[3]-5*LinearModel.2$coef[4] (Intercept) 79.5399 > LinearModel.2$coef[1]+20*LinearModel.2$coef[2]+ 4*LinearModel.2$coef[3]+5*LinearModel.2$coef[4] (Intercept) 95.83004 Les autres variables restant identiques, il semble donc que la concentration dozone augmente lorsque la valeur de la variable Vx augmente. On sy attendait puisque le coecient 3 est positif. Remarque 5 Cette question permet dinsister sur un point important : en rgression multiple, on e ninterpr`te la valeur dun coecient associ a une variable quen considrant les autres variables e e ` e xes. Ainsi, si un coecient j particulier est positif, cela ne veut pas dire ncssairement quune e e e augmentation de la variable Xj produit une augmentation de Y . On peut par contre dire quelle produit une augmentation de Y , les autres variables tant par ailleurs xes. e e

LinearModel.2$residuals

40 0

20

20

40

20

40

60 Index

80

100

Fig. 8 Reprsentations des rsidus. e e

Exercice 8 : Choix de mod`le e

1. Le jeu de donnes contient 10 variables explicatives quantitatives : T9, T12, T15, Ne9, Ne12, Ne15, e Vx9, Vx12, Vx15 et maxO3v. 2. Le mod`le complet scrit donc e e maxO3i =0 + 1 T 9i + 2 T 12i + 3 T 15i + 4 N e9i + 5 N e12i + 6 N e15i + 7 V x9i + 8 V x12i + 9 V x15i + 10 maxO3vi + i . 3. Le mod`le complet ajust est le suivant : e e > LinearModel.1 <- lm(maxO3 ~ T9 +T12 +T15 +Ne9 +Ne12 +Ne15 +Vx9 + Vx12 +Vx15 + +maxO3v, data=Dataset) > summary(LinearModel.1) Call: lm(formula = maxO3 ~ T9 + T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx12 + Vx15 + maxO3v, data = Dataset) Residuals: Min 1Q -53.566 -8.727

Median -0.403

3Q 7.599

Max 39.458

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 12.24442 13.47190 0.909 0.3656

T9 -0.01901 T12 2.22115 T15 0.55853 Ne9 -2.18909 Ne12 -0.42102 Ne15 0.18373 Vx9 0.94791 Vx12 0.03120 Vx15 0.41859 maxO3v 0.35198 --Signif. codes: 0 ***

1.12515 1.43294 1.14464 0.93824 1.36766 1.00279 0.91228 1.05523 0.91568 0.06289

-0.017 0.9866 1.550 0.1243 0.488 0.6266 -2.333 0.0216 * -0.308 0.7588 0.183 0.8550 1.039 0.3013 0.030 0.9765 0.457 0.6486 5.597 1.88e-07 ***

0.001 ** 0.01 * 0.05 . 0.1 1

Residual standard error: 14.36 on 101 degrees of freedom Multiple R-squared: 0.7638,Adjusted R-squared: 0.7405 F-statistic: 32.67 on 10 and 101 DF, p-value: < 2.2e-16 4. Au seuil 5% seuls les coecients associs aux variables maxO3v et Ne9 sont signicatifs. e 5. Le mod`le linaire ajust avec ces deux variables donne : e e e > LinearModel.2 <- lm(maxO3 ~ Ne9 + maxO3v, data=Dataset) > summary(LinearModel.2) Call: lm(formula = maxO3 ~ Ne9 + maxO3v, data = Dataset) Residuals: Min 1Q -62.484 -9.631

Median -2.208

3Q 10.315

Max 53.077

Coefficients: Estimate Std. Error t value (Intercept) 65.25915 6.86441 9.507 Ne9 -5.08592 0.62098 -8.190 maxO3v 0.55327 0.05699 9.709 --Signif. codes: 0 *** 0.001 ** 0.01

Pr(>|t|) 5.71e-16 *** 5.43e-13 *** < 2e-16 *** * 0.05 . 0.1 1

Residual standard error: 16.32 on 109 degrees of freedom Multiple R-squared: 0.671,Adjusted R-squared: 0.665 F-statistic: 111.2 on 2 and 109 DF, p-value: < 2.2e-16 6. Si on rajoute la variable T12, on obtient : > LinearModel.3 <- lm(maxO3 ~ Ne9 + maxO3v +T12, data=Dataset) > summary(LinearModel.3) Call: lm(formula = maxO3 ~ Ne9 + maxO3v + T12, data = Dataset) Residuals: Min

1Q

Median

3Q

Max

-56.385

-7.872

-1.941

7.899

41.513

Coefficients: Estimate Std. Error t value (Intercept) 9.76225 11.10038 0.879 Ne9 -3.02423 0.64342 -4.700 maxO3v 0.37571 0.05801 6.477 T12 2.85308 0.48052 5.937 --Signif. codes: 0 *** 0.001 ** 0.01

Pr(>|t|) 0.381 7.71e-06 *** 2.85e-09 *** 3.57e-08 *** * 0.05 . 0.1 1

Residual standard error: 14.23 on 108 degrees of freedom Multiple R-squared: 0.752,Adjusted R-squared: 0.7451 F-statistic: 109.1 on 3 and 108 DF, p-value: < 2.2e-16 On remarque que le coecient associ ` cette nouvelle variable est signicatif. On comprend que ne ea retenir que les variables associes aux coecents signicatifs du mod`le complet ne donne pas une e e procdure de choix de mod`le satisfaisante. e e 7. Les nuages de points croisant toutes les variables explicatives sobtient ainsi (cf. gure 9) : > [Link](~maxO3v+Ne9+Ne12+Ne15+T9+T12+T15+Vx9+Vx12+Vx15, + [Link]=FALSE, smooth=FALSE, span=0.5, diagonal = density, data=Dataset) On voit que certaines variables sont tr`s corrles entre elles comme les mesures de temprature. e ee e Les maintenir toutes les trois produit une dilution de linformation qui a tendance ` masquer leur a eet. Sans doute faut-il ne retenir que lune des trois par exemple. Procdure pas ` pas On enl`ve successivement les variables : T9, Vx12, Ne15, Ne12, Vx15 et T15 e a e pour retenir le mod`le ` 4 variables : e a > summary(LinearModel.1) Call: lm(formula = maxO3 ~ T12 + Ne9 + Vx9 + maxO3v, data = Dataset) Residuals: Min 1Q -52.396 -8.377

Median -1.086

3Q 7.951

Max 40.933

Coefficients: Estimate Std. Error t value (Intercept) 12.63131 11.00088 1.148 T12 2.76409 0.47450 5.825 Ne9 -2.51540 0.67585 -3.722 Vx9 1.29286 0.60218 2.147 maxO3v 0.35483 0.05789 6.130 --Signif. codes: 0 *** 0.001 ** 0.01

Pr(>|t|) 0.253443 6.07e-08 0.000317 0.034055 1.50e-08

*** *** * ***

* 0.05 . 0.1 1

Residual standard error: 14 on 107 degrees of freedom Multiple R-squared: 0.7622,Adjusted R-squared: 0.7533 F-statistic: 85.75 on 4 and 107 DF, p-value: < 2.2e-16

15

8 2 160

maxO3v
||||||||||||||||||| ||||||| || |||| || || |||| |||| || |

Ne9
| | | | | || | | || || | | || || ||

Ne12
|| | | || | | | | | || | ||| || | ||

Ne15
||||||||| || || || || || |||

T9
|||||||||||| ||| ||||| || || |||| ||| ||| |||

T12
15 |||||||||||||| || |||||||||||| |||| ||||| ||||||| |||||| || ||| |

T15
|||||||||||||||| |||||||||||| |||||||| || |||||||| ||| ||| ||

Vx9
|||||||||||||||||| |||| |||| || ||| | || 5 5

Vx12
||||||||||||||||| ||||||||||| ||||||| | |||||| ||| |||| 8 2

Vx15
||||||||||||||||| |||||||| |||||||| || || |||| 40 160 0 6 15 15 35 5 5

Fig. 9 Nuages des points en variables.

Exercice 9 : Autre exemple

On peut se placer sous Rcmdr pour importer le jeu de donnes. Il faut ici spcier que le sparateur de e e e champ est : et la dcimale est dnie par ,. e e > Dataset <+ [Link]("/home/jegou/ENSEIGNEMENT/GEO/MASTER/COURS-TD/S2/[Link]", + header=TRUE, sep=":", [Link]="NA", dec=",", [Link]=TRUE) On part du mod`le complet ` 12 variables explicatives. On enl`ve la variable associe au coecient le moins e a e e signicatif ; on recalcule les coecients associs ` ce mod`le puis on enl`ve ` nouveau la variable la moins e a e e a signicative... jusqu` ce que tous les coecients soient signicatifs (au seuil 5% par exemple). a A lissue de cette procdure, on conserve le mod`le suivant impliquant 4 variables : e e > LinearModel.1 <- lm(Y ~ X1 +X2 +X5 +X6 ,

15

35

15

40

data=Dataset)

> summary(LinearModel.1) Call: lm(formula = Y ~ X1 + X2 + X5 + X6, data = Dataset) Residuals: Min 1Q -53.122 -13.129

Median 1.264

3Q 13.964

Max 57.050

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -14.1011 2.9437 -4.79 2.64e-06 *** X1 2.9863 0.2220 13.45 < 2e-16 *** X2 9.3518 0.2179 42.92 < 2e-16 *** X5 -6.9478 0.2076 -33.46 < 2e-16 *** X6 12.4763 0.2316 53.87 < 2e-16 *** --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 20.41 on 295 degrees of freedom Multiple R-squared: 0.9541,Adjusted R-squared: 0.9535 F-statistic: 1533 on 4 and 295 DF, p-value: < 2.2e-16

Analyse de covariance
1. Par linstruction suivante, on obtient la gure 10 : > scatterplot(maxO3~T12 | pluie, [Link]=lm, smooth=FALSE, labels=FALSE, + boxplots=FALSE, span=0.5, [Link]=TRUE, data=Dataset) 2. On obtient les param`tres de la droite de rgression des donnes par temps sec de la mani`re suivante : e e e e > RegModel.1 <- lm(maxO3~T12, data=Dataset, subset=pluie=="Sec") > summary(RegModel.1) Call: lm(formula = maxO3 ~ T12, data = Dataset, subset = pluie == "Sec") Residuals: Min 1Q -36.3567 -11.9778

Exercice 10 : Ecriture de tous les mod`les e

Median 0.1406

3Q 11.5963

Max 42.4064

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -21.2601 12.2970 -1.729 0.0884 . T12 5.3255 0.5278 10.091 4.44e-15 *** --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

pluie Pluie Sec 160 maxO3 40 60 80 100 120 140

15

20 T12

25

30

Fig. 10 Droites de rgression dans chaque groupe de pluie. e Residual standard error: 18.19 on 67 degrees of freedom Multiple R-squared: 0.6031,Adjusted R-squared: 0.5972 F-statistic: 101.8 on 1 and 67 DF, p-value: 4.445e-15 Par temps de pluie : > RegModel.2 <- lm(maxO3~T12, data=Dataset, subset=pluie=="Pluie") > summary(RegModel.2) Call: lm(formula = maxO3 ~ T12, data = Dataset, subset = pluie == "Pluie") Residuals: Min 1Q -27.137 -10.607

Median -1.630

3Q 7.565

Max 43.319

Coefficients: Estimate Std. Error t value (Intercept) 7.0386 17.5636 0.401 T12 3.4419 0.9033 3.810 --Signif. codes: 0 *** 0.001 ** 0.01

Pr(>|t|) 0.690687 0.000457 *** * 0.05 . 0.1 1

Residual standard error: 14.94 on 41 degrees of freedom Multiple R-squared: 0.2615,Adjusted R-squared: 0.2435 F-statistic: 14.52 on 1 and 41 DF, p-value: 0.0004574

3. On observe une pente plus forte lorsque le temps est sec : on peut donc penser que que les variations de maxO3 sont plus sensibles aux augmentations de temprature par temps sec que par temps pluvieux. e 4. Le mod`le complet scrit donc : e e maxO3i = 0,p + 1,p T 12i + i maxO3i = 0,s + 1,s T 12i + i On lcrit sur R de la mani`re suivante : e e > LinearModel.3 <- lm(maxO3 ~ T12:pluie +pluie - 1, data=Dataset) > summary(LinearModel.3) Call: lm(formula = maxO3 ~ T12:pluie + pluie - 1, data = Dataset) Residuals: Min 1Q -36.35669 -11.04905 Coefficients: Estimate Std. Error t value Pr(>|t|) pluiePluie 7.0386 20.0222 0.352 0.72587 pluieSec -21.2601 11.5118 -1.847 0.06751 . T12:pluiePluie 3.4419 1.0298 3.342 0.00114 ** T12:pluieSec 5.3255 0.4941 10.779 < 2e-16 *** --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 17.03 on 108 degrees of freedom Multiple R-squared: 0.9687,Adjusted R-squared: 0.9676 F-statistic: 836.6 on 4 and 108 DF, p-value: < 2.2e-16 Remarque 6 Calculer les param`tres du mod`le complet revient a calculer les quations des deux e e ` e droites de rgression prcdentes. e e e 5. Le mod`le avec une pente commune mais des ordonnes ` lorigine direntes scrit : e e a e e maxO3i = 0,p + 1 T 12i + i maxO3i = 0,s + 1 T 12i + i Sur R, cela donne : > LinearModel.4 <- lm(maxO3 ~ T12 +pluie - 1, data=Dataset) > summary(LinearModel.4) Call: lm(formula = maxO3 ~ T12 + pluie - 1, data = Dataset) Residuals: Min 1Q -36.2009 -12.4102

Median 0.05153

3Q 10.72335

Max 43.31871

Median -0.3394

3Q 9.5577

Max 44.9709

Coefficients: Estimate Std. Error t value Pr(>|t|) T12 4.973 0.449 11.077 <2e-16 *** pluiePluie -22.480 9.042 -2.486 0.0144 * pluieSec -13.179 10.499 -1.255 0.2121 --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 17.16 on 109 degrees of freedom Multiple R-squared: 0.9679,Adjusted R-squared: 0.9671 F-statistic: 1097 on 3 and 109 DF, p-value: < 2.2e-16 6. Le mod`le avec une seule ordonne ` lorigine mais des pentes direntes scrit lui : e e a e e maxO3i = 0 + 1,p T 12i + i maxO3i = 0 + 1,s T 12i + i Sur R : > LinearModel.5 <- lm(maxO3 ~ T12:pluie, data=Dataset) > summary(LinearModel.5) Call: lm(formula = maxO3 ~ T12:pluie, data = Dataset) Residuals: Min 1Q -36.43933 -11.63807 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -14.2295 10.0028 -1.423 0.158 T12:pluiePluie 4.5265 0.5274 8.583 7.13e-14 *** T12:pluieSec 5.0286 0.4316 11.652 < 2e-16 *** --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 17.07 on 109 degrees of freedom Multiple R-squared: 0.64,Adjusted R-squared: 0.6334 F-statistic: 96.89 on 2 and 109 DF, p-value: < 2.2e-16 7. Le mod`le qui ne prend pas en compte la variable pluie scrit : e e maxO3i = 0 + 1 T 12i + i Remarque 7 Cela revient donc a ajuster la droite de rgression toutes donnes de la variable pluie ` e e confondues comme en rgression simple. e On retrouve donc les rsultats de lexercice dapplication en section 1 : e > LinearModel.6 <- lm(maxO3 ~ T12, data=Dataset) > summary(LinearModel.6) Call:

Median -0.06564

3Q 9.82501

Max 44.84680

lm(formula = maxO3 ~ T12, data = Dataset) Residuals: Min 1Q -38.0789 -12.7352

Median 0.2567

3Q 11.0029

Max 44.6714

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -27.4196 9.0335 -3.035 0.003 ** T12 5.4687 0.4125 13.258 <2e-16 *** --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 17.57 on 110 degrees of freedom Multiple R-squared: 0.6151,Adjusted R-squared: 0.6116 F-statistic: 175.8 on 1 and 110 DF, p-value: < 2.2e-16 8. Enn, le mod`le qui ne tient pas compte de la variable T12 scrit : e e maxO3i = 0,p + i maxO3i = 0,s + i Remarque 8 On retrouve le mod`le danalyse de variance a un facteur. e ` Sur R, cela donne : > LinearModel.7 <- lm(maxO3 ~ pluie, data=Dataset) > summary(LinearModel.7) Call: lm(formula = maxO3 ~ pluie, data = Dataset) Residuals: Min 1Q -37.841 -17.841 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 73.395 3.798 19.324 < 2e-16 *** pluie[[Link]] 27.445 4.839 5.672 1.16e-07 *** --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 24.91 on 110 degrees of freedom Multiple R-squared: 0.2263,Adjusted R-squared: 0.2192 F-statistic: 32.17 on 1 and 110 DF, p-value: 1.157e-07

Median -4.395

3Q 11.409

Max 65.159

Exercice 11 : Choix du meilleur mod`le e

Les mod`les dnis sont les suivants : e e # modele complet : > LinearModel.3 <- lm(maxO3 ~ T12:pluie +pluie - 1, data=Dataset)

# modele avec pentes egales LinearModel.4 <- lm(maxO3 ~ T12 +pluie - 1, data=Dataset) # modele avec ordonnees a lorigine egales LinearModel.5 <- lm(maxO3 ~ T12:pluie, data=Dataset) # modele sans effet pluie LinearModel.6 <- lm(maxO3 ~ T12, data=Dataset) # modele danalyse de variance LinearModel.7 <- lm(maxO3 ~ pluie, data=Dataset) Nous testons le mod`le complet (3) contre le mod`le avec pentes gales (4) : e e e > anova(LinearModel.3, LinearModel.4) Analysis of Variance Table Model 1: Model 2: [Link] 1 108 2 109 maxO3 ~ T12:pluie + pluie - 1 maxO3 ~ T12 + pluie - 1 RSS Df Sum of Sq F Pr(>F) 31313 32102 -1 -789 2.7197 0.1020

Nous retenons le mod`le (4) avec pentes gales. e e Testons maintenant le mod`le (4) contre le simple mod`le danalyse de variance (7) : e e > anova(LinearModel.4, LinearModel.7) Analysis of Variance Table Model 1: maxO3 ~ T12 + pluie - 1 Model 2: maxO3 ~ pluie [Link] RSS Df Sum of Sq F Pr(>F) 1 109 32102 2 110 68238 -1 -36136 122.70 < 2.2e-16 *** La probabilit critique est faible, il vaut mieux le mod`le avec pentes gales et ordonnes ` lorigine direntes e e e e a e que le simple mod`le danalyse de variance (qui reviendrait ` considrer que la pente commune est nulle). e a e Testons le mod`le (4) contre le mod`le de rgression simple (6) : e e e > anova(LinearModel.4, LinearModel.6) Analysis of Variance Table Model 1: maxO3 ~ T12 + pluie - 1 Model 2: maxO3 ~ T12 [Link] RSS Df Sum of Sq F Pr(>F) 1 109 32102 2 110 33948 -1 -1846 6.2689 0.01377 * --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Le mod`le (4) avec deux ordonnes ` lorigine est meilleur. e e a Cest donc le mod`le (4) que nous retenons. e