MAT 3775 - Analyse de la régression

• Télécharger R de CRAN, voir http://www.r-project.org/.
  • Suivre le lien pour télécharger R.
  • Choisir un des sites miroir, par exemple Toronto ou Montréal.
  • Sélectionner votre système d'exploitation, par exemple Windows, Mac OS X ou Linux.
  • Télécharger et installer la base.
• Voici un bon vidéo de You Tube qui illustre le téléchargement et l'installation de R : https://www.youtube.com/watch?v=jM1KDt5gMVM.

Exemple 1 : Travailler à la console

Ci-dessous, on a un affiche de la console de R. Nous entrons des commandes à l'invite >.

Nous pouvons utiliser R comme une calculatrice. Nous avons calculé (1/8)*ln(2^3)+5. La réponse est 5,25993.

Il y a une console dans Rstudio. C'est dans la fenêtre en bas à gauche.

x=c(16,15,14,16,13,12,14,13,10)

y<-c(16,15,14,16,13,12,14,13,10)
x<-c(0,0,0,25,25,25,50,50,50)

> 1-pnorm(1.65)
[1] 0.04947147
> 1-pchisq(12.6,6)
[1] 0.04984649
> 1-pt(1.71,24)
[1] 0.05008269
> 1-pf(1.85,20,15)
[1] 0.1140452

> qt(0.975,43)
[1] 2.016692

> rnorm(10,46,7)
 [1] 46.93336 55.43874 39.04338 56.84523 52.93618 39.41287 31.94741 43.18834
 [9] 45.07614 51.84981

• Considérons une étude avec des volontaires en bonne santé. Un stimulus est appliqué aux doigts du sujet et on mesure la vitesse de conduction de la moelle épinière (VC). On veut décrire l'association entre la taille de l'indivu (en cm) et la vitesse de conduction de la moelle épinière pour les individus en bonne santé. On a utilise les commandes suivantes pour importer les données et produire le nuage de points de la taille contre la vitesse de condution de la moelle épinière.
• Nous allons sauvegarder les données dans notre répertoire de travail (current working directory). Pour afficher le répertoire de travail, on utilise getwd(). On peut modifier le répertoire avec la fonction setwd().
  
  Voici des examples de l'utilisation de la fonction setwd().
  
  

  setwd(directory)      # change to directory
  setwd("c:/Rstuff/directory")  # note / instead of \ in windows
  setwd("/usr/gilles/directory") # on linux or a Mac
  

• Les données sont dans le fichier suivant: VC.csv. Sauvegarder le fichier dans votre répertoire de travail (current working directory).
  
  On peut importer un jeu des données d'un fichier texte avec la commande suivante.

  VC <- read.csv("VC.csv")
  

  Remarques:
  • La fonction ci-haut a produit un objet qui est un dataframe (c-à-d un jeu de données). Le nom de l'objet est VC.
  • Un jeu de données (c-à-d un dataframe) est un tableau de données. Une rangée est une unité statistique. Les variables (c-à-d les colonnes) sont utilisées pour décrire les unités.
  On affiche les noms des colonnes avec la fonction names.

  > names(VC)
  [1] "Taille.en.cm" "VC"
  

  On peut référer à une colonne avec son nom ou son indice.

  > VC$Taille.en.cm
    [1] 149 149 155 155 156 156 157 157 158 158 160 160 161 161 161 161 161 161
   [19] 162 162 168 168 168 168 168 168 168 168 170 170 170 170 171 171 171 171
   [37] 172 172 172 172 181 181 181 182 182 182 182 182 182 182 184 184 184 184
   [55] 184 184 185 185 187 187 162 162 163 163 163 163 163 163 163 164 164 164
   [73] 164 165 165 165 165 165 166 166 166 166 166 167 167 167 167 167 167 167
   [91] 167 168 173 173 174 174 175 175 175 175 175 175 175 175 176 176 176 176
  [109] 176 177 177 179 179 179 179 179 179 179 179 179 180 180 181 181 187 187
  [127] 188 188 189 189 190 190 190 190 191 191 191 191 194 194 194 194 195 195
  [145] 196 196 197 197 200 200 202 202 182 190 190
  > VC[,1]
    [1] 149 149 155 155 156 156 157 157 158 158 160 160 161 161 161 161 161 161
   [19] 162 162 168 168 168 168 168 168 168 168 170 170 170 170 171 171 171 171
   [37] 172 172 172 172 181 181 181 182 182 182 182 182 182 182 184 184 184 184
   [55] 184 184 185 185 187 187 162 162 163 163 163 163 163 163 163 164 164 164
   [73] 164 165 165 165 165 165 166 166 166 166 166 167 167 167 167 167 167 167
   [91] 167 168 173 173 174 174 175 175 175 175 175 175 175 175 176 176 176 176
  [109] 176 177 177 179 179 179 179 179 179 179 179 179 180 180 181 181 187 187
  [127] 188 188 189 189 190 190 190 190 191 191 191 191 194 194 194 194 195 195
  [145] 196 196 197 197 200 200 202 202 182 190 190
  

  
  
  Commentaire :
  • [,1] veut dire toutes les rangés, mais seulement la colonne 1.

> cerveau<-read.csv("cerveaux.csv")
> head(cerveau)
  masse..grammes. gestation..en.jours.
1              10                   50
2              10                   55
3              12                   56
4              12                   75
5              12                   85
6              45                   86

with(cerveau, plot(gestation..en.jours.,masse..grammes.,xlab="Gestation (en jours)",ylab="Masse du cerveau (en grammes)"))

## ajuster une courbe lisse (loess=lowess=locally weighted scatterplot smoothing)
mod.loess<-loess(masse..grammes.~gestation..en.jours.,data=cerveau)
## obtenir l'étendue pour x
range(cerveau$gestation..en.jours.)
## construire un nouveau jeu de données
new.data<-data.frame(gestation..en.jours.=seq(50,270,by=10))
## add Lowess Smooth to the plot
with(new.data,lines(x=gestation..en.jours.,y=predict(mod.loess,new.data),lty=2))

• La commande lines(x,y) ajoute au nuage de points un nuage de points de y contre x tel qu'on a des segments de droite entre le points successifs.
• L'argument lty décrit le type de ligne qui sera utilisé.
• Avec la fonction data.frame, on a construit un nouveau jeu de données qui contient une colonne qui a le même nom que notre variable x.
• La fonction predict(object, newdata) évalue un modèle estimé object aux valeurs dans le jeu de données newdata.

• La corrélation entre la masse du cerveau et la durée de la période de gestation est $r=0.81$. Ceci est la corrélation de Pearson. Elle est une description de l'intensité d'une association linéaire entre deux variables numériques. Ici le lien n'est pas linéaire, alors c'est difficile à interpréter la corrélation. C'est mieux d'utiliser la corrélation de Spearmen quand le lien est monotone et non-linéaire.
  
  

  > with(cerveau,cor(gestation..en.jours.,masse..grammes.))
  [1] 0.8113314
  

  
  
  
• Par défaut, la fonction method="pearson". On peut plutôt utiliser l'argument method="spearman" pour obtenir la corrélation de Spearman. (C'est la corrélation de Pearson entre les rangs de x et les rangs de y.)
• La corrélation entre la masse du cerveau et la durée de la période de gestation est $r_S=0.86$.
  
  

  > with(cerveau,cor(gestation..en.jours.,masse..grammes.,method="spearman"))
  [1] 0.8628308
  

# estimer la droite
mod<-lm(masse..grammes.~gestation..en.jours.,data=cerveau)
# ajouter la droite au diagramme
abline(mod)
legend("topleft",legend=c("droite estimée","courbe lisse"),
       lty=c(1,2),inset=0.05)

a+x        # additionne a à chaque composante
a*x        # multiplier a à chaque composante
x/a        # diviser chaque composante par a
x+y        # addition par composante
x*y        # multiplication par composante
log(x)     # appliquer ln à chaque composante
sqrt(x)    # appliquer la racine carrée à chaque composante
x^a        # appliquer la puissance à chaque composante
abs(x)     # appliquer la valeur absolue à chaque composante

sum(x)       # la somme des composantes de x
length(x)    # le nombre de composantes

> x=c(10, 12, 11, 12, 15, 20, 16)
> sum((x-mean(x))^2)/(length(x)-1)
[1] 12.2381
> var(x)
[1] 12.2381

> sum((x-mean(x))*(y-mean(y)))/(length(x)-1)
[1] 9.404762
> cov(x,y)
[1] 9.404762

> carburant<-read.csv("evaporation.csv")
> head(carburant)
  Velocite Evaporation
1       20        0.18
2       60        0.37
3      100        0.35
4      140        0.78
5      180        0.56
6      220        0.75

> x<-carburant$Velocite
> y<-carburant$Evaporation
> sum(x)
[1] 2220
> sum(y)
[1] 9.1
> sum(x^2)
[1] 580400
> sum(y^2)
[1] 9.6722
> sum(x*y)
[1] 2340.4

> carburant<-read.csv("evaporation.csv")
> head(carburant)
  Velocite Evaporation
1       20        0.18
2       60        0.37
3      100        0.35
4      140        0.78
5      180        0.56
6      220        0.75

> mod<-lm(Evaporation ~ Velocite, data=carburant)
> mod

Call:
lm(formula = Evaporation ~ Velocite, data = carburant)

Coefficients:
(Intercept)     Velocite
   0.059033     0.003807

• Nous avons construit un objet lm que nous avons nommé mod. Cet objet contient a plusieurs composantes qu'on peut afficher en utilisant la fonction names.

  > names(mod)
   [1] "coefficients"  "residuals"     "effects"       "rank"          "fitted.values" "assign"        "qr"
   [8] "df.residual"   "xlevels"       "call"          "terms"         "model"
  

  
  
  La première composantes est le vecteur des esimations des betas. Alors, $b_0=0.059032967$ et $b_1= 0.003806593$.

  > mod$coefficients
  (Intercept)    Velocite
  0.059032967 0.003806593
  

  
  
  La deuxième composante est le vecteur des résidus et la 8ième composante est le nombre de degrés de liberté des résidus. On peut les utilisés pour calculer l'estimation de la variance de l'erreur, c'est-à-dire $MSE=\sum_{i=1}^n e_i^2/(n-2) =0.02512$. L'écart type de cette variance est dite l'écart type résiduelle $s_e=\sqrt{\mbox{MSE}}=0.1584~\mbox{mm}^2/s$. Ceci est une description d'un écart type de la droite.

  > MSE<-sum(mod$residuals^2)/mod$df.residual
  > MSE; sqrt(MSE)
  [1] 0.02511653
  [1] 0.158482
  

  
  
  
• Il y a plusieurs fonctions qu'on puisse utiliser avec un object lm. Pour afficher ces fonctions, on peut utiliser la fonction methods.

  > methods(class=lm)
   [1] add1           alias          anova          case.names     coerce
   [6] confint        cooks.distance deviance       dfbeta         dfbetas
  [11] drop1          dummy.coef     effects        extractAIC     family
  [16] formula        hatvalues      influence      initialize     kappa
  [21] labels         logLik         model.frame    model.matrix   nobs
  [26] plot           predict        print          proj           qr
  [31] residuals      rstandard      rstudent       show           simulate
  [36] slotsFromS3    summary        variable.names vcov
  see '?methods' for accessing help and source code
  

  
  
  Voici une description de l'ajustment avec le maximum du log la fonction de vraisemblance $$\ell=-\frac{n}{2}\ln(2\,\pi\,\mbox{SSE}/n)-n/2=5,\!758624.$$ Il y a 3 degrés de liberté, puisque nous avons 3 paramètres à estimer $\beta_0,\beta_1$ et $\sigma^2$.

  > logLik(mod)
  'log Lik.' 5.758624 (df=3)
  

  
  
  V\'erifions que R ait bien utilis\'e la formule ci-haut pour calculer le maximum du log la fonction de vraisemblance.

  > ## p=2=le nombre de paramètres dans la fonction de la moyenne
  > p<-length(mod$coefficients)
  > ## n = (n-p) + p = degrés de liberté de l'erreur + p
  > n<-mod$df.residual+p
  > ## SSE = somme des carrés résiduelle
  > SSE<-sum(mod$residuals^2)
  > # calculons le max du log de la vraisemblance
  > -n/2*log(2*pi*SSE/n)-n/2
  [1] 5.758624
  

• On peut obtenir des intervalles de confiance pour les estimations des param\`etres.

  > confint(mod)
                    2.5 %      97.5 %
  (Intercept) -0.16731971 0.285385640
  Velocite     0.00282118 0.004792006
  

  
  
  Commentaires:
  • Le niveau de confiance est 95% par défaut. Mais, on peut modifier l'argument level pour changer le niveau de confiance.

    > confint(mod,level=0.98)
                         1 %        99 %
    (Intercept) -0.223281643 0.341347577
    Velocite     0.002577553 0.005035633
    

  • On peut aussi spécifier les paramètres qui nous intéressent. Supposons qu'on veut seulement un IC pour la pente.

    > confint(mod,parm=c("Velocite"))
                  2.5 %      97.5 %
    Velocite 0.00282118 0.004792006
    

• On peut obtenir un sommaire de l'ajustement avec la fonction summary.

  > summary(mod)
  
  Call:
  lm(formula = Evaporation ~ Velocite, data = carburant)
  
  Residuals:
       Min       1Q   Median       3Q      Max
  -0.18422 -0.14648  0.04483  0.13786  0.18804
  
  Coefficients:
               Estimate Std. Error t value Pr(>|t|)
  (Intercept) 0.0590330  0.1000605   0.590     0.57
  Velocite    0.0038066  0.0004356   8.739 1.09e-05 ***
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  Residual standard error: 0.1585 on 9 degrees of freedom
  Multiple R-squared:  0.8946,    Adjusted R-squared:  0.8829
  F-statistic: 76.36 on 1 and 9 DF,  p-value: 1.086e-05
  

  
  
  Commentaire:
  • Ce n'est pas nécessaire d'afficher toute le sommaire. Le sommaire lui-même est un objet avec des composantes.

    > names(summary(mod))
     [1] "call"          "terms"         "residuals"     "coefficients"  "aliased"       "sigma"
     [7] "df"            "r.squared"     "adj.r.squared" "fstatistic"    "cov.unscaled"
    

    Voici l'estimations des paramètres et le test de la signification pour chacun.

    > summary(mod)$coefficients
                   Estimate   Std. Error   t value     Pr(>|t|)
    (Intercept) 0.059032967 0.1000605427 0.5899725 5.697228e-01
    Velocite    0.003806593 0.0004356076 8.7385826 1.085924e-05
    

    On peut afficher le coefficient de détermination $R^2=89,\!5\%$.

    > summary(mod)$r.squared
    [1] 0.8945677
    

    On peut afficher l'écart type résiduelle $s_e=\sqrt{\mbox{MSE}}=0,\!158482$.

    > summary(mod)$sigma
    [1] 0.158482
    

> VC <- read.csv("VC.csv")
> head(VC)
  Taille.en.cm   VC
1          149 14.4
2          149 13.4
3          155 13.5
4          155 13.5
5          156 13.0
6          156 13.6

> mod<-lm(VC ~ Taille.en.cm, data=VC)
> anova(mod)
Analysis of Variance Table

Response: VC
              Df Sum Sq Mean Sq F value    Pr(>F)
Taille.en.cm   1 287.56 287.563   391.3 < 2.2e-16 ***
Residuals    153 112.44   0.735
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

• La statistique du test est $F=391,\!3$. Ceci veut dire que l'estimation de la variance de l'erreur basée sur la somme des carrés de la régression est 391,\!3 fois aussi grande que l'estimation de la variance de l'erreur basée sur la somme des carrés résiduelle. Il est fort probable que $\mbox{MSR}$ n'estime pas $\sigma^2$, mais plut\^ot quelque chose de plus grand. On peut d\'emontrer que $$E\{\mbox{MSR}\}=E\{\mbox{SSR}\}=\sigma^2\,(1+\beta_1^2\,s_{xx}).$$ Ceci veut dire que nous avons des preuves que $\beta_1\neq 0$.
• Pour avoir une mesure de la signification des preuves contre $H_0:\beta_1=0$, on calcul nos chances d'avoir observer une statistique $F$ aussi grande que $391,\!3$ en supposant que $H_0$ est vraie. Sous $H_0$, $F^*\sim F(1,n_2)$, alors la valeur $P$ est $$P(F(1,153)> 391,\!3) < 0,\!0001.$$
• Le test basé sur $T$ et l'analyse de variance sont équivalents, puisqu'on peut démontrer que $$t^*=\left( \frac{b_1}{s\{b_1\}} \right)^2=\frac{\mbox{SSR}/1}{\mbox{SSE}/(n-2)}=F ~~~~ \mbox{ et } ~~~~ t_{1-\alpha/2;n-2}^2=F(1-alpha;1,n-2).$$

> VC <- read.csv("VC.csv")
> head(VC)
  Taille.en.cm   VC
1          149 14.4
2          149 13.4
3          155 13.5
4          155 13.5
5          156 13.0
6          156 13.6

> mod<-lm(VC ~ Taille.en.cm, data=VC)
> ## modele réduit
> mod0<-lm(VC ~ 1, data=VC)
> mod<-lm(VC ~ Taille.en.cm, data=VC)
> ## modele réduit
> mod0<-lm(VC ~ 1, data=VC)
> anova(mod0,mod)
Analysis of Variance Table

Model 1: VC ~ 1
Model 2: VC ~ Taille.en.cm
  Res.Df    RSS Df Sum of Sq     F    Pr(>F)
1    154 400.00
2    153 112.44  1    287.56 391.3 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> azote<-read.csv("Azote.csv")
> head(azote)
  Masse Azote
1  0.59   non
2  0.47   non
3  0.25   non
4  0.36   non
5  0.42   non
6  0.19   non

> levels(azote$Azote)
[1] "non" "oui"
> azote$Azote<-factor(azote$Azote, levels=c("oui","non"))
> levels(azote$Azote)
[1] "oui" "non"

> azote$Azote<-factor(azote$Azote, levels=c("non","oui"))

> aggregate(Masse~Azote,data=azote,FUN=mean)
  Azote Masse
1   non 0.392
2   oui 0.628
> aggregate(Masse~Azote,data=azote,FUN=var)
  Azote      Masse
1   non 0.01355111
2   oui 0.01819556
> aggregate(Masse~Azote,data=azote,FUN=length)
  Azote Masse
1   non    10
2   oui    10

> moyenne<-with(azote,tapply(Masse,Azote,mean))
> variance<-with(azote,tapply(Masse,Azote,var))
> n<-with(azote,tapply(Masse,Azote,length))
> data.frame(moyenne,variance,n)
    moyenne   variance  n
non   0.392 0.01355111 10
oui   0.628 0.01819556 10

> ## variance poindérée
> sum(variance*(n-1)/sum(n))
[1] 0.014286

boxplot(Masse~Azote,data=azote,ylab="Masse",
names=c("Sans azote","Avec azote"))
stripchart(Masse~Azote,data=azote,method="jitter",
vertical = TRUE,add=TRUE)

> t.test(Masse~Azote,data=azote,var.equal=TRUE)

        Two Sample t-test

data:  Masse by Azote
t = -4.1885, df = 18, p-value = 0.0005521
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.3543747 -0.1176253
sample estimates:
mean in group non mean in group oui
            0.392             0.628

> contrasts(azote$Azote)
    oui
non   0
oui   1

> mod<-lm(Masse~Azote,data=azote)
> mod$coefficients
(Intercept)    Azoteoui
      0.392       0.236
> summary(mod)$sigma^2
[1] 0.01587333

> summary(mod)$coefficients
            Estimate Std. Error  t value     Pr(>|t|)
(Intercept)    0.392 0.03984135 9.839024 1.145796e-08
Azoteoui       0.236 0.05634418 4.188543 5.520945e-04

> VC <- read.csv("VC.csv")
> head(VC)
  Taille.en.cm   VC
1          149 14.4
2          149 13.4
3          155 13.5
4          155 13.5
5          156 13.0
6          156 13.6

> with(VC,cor.test(Taille.en.cm,VC))

        Pearson's product-moment correlation

data:  Taille.en.cm and VC
t = 19.781, df = 153, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.7967316 0.8869721
sample estimates:
      cor
0.8478829

• Pour obtenir l'intervalle de confiance pour la corrélation, la fonction cor.test utilise la transformation $z$ de Fisher. Cette transformation est $$\frac{1}{2}\ln\left( \frac{1+r}{1-r} \right)=\mbox{arctanh}(r).$$ Fisher a démontré que la distribution d'échantillonage de cette transformation suit approximativement $$N\left(\frac{1}{2}\ln\left( \frac{1+\rho}{1-\rho} \right), \frac{1}{n-3} \right).$$ Alors, un IC \`a 95\% pour $\mbox{arctanh}(\rho)$ est $$\frac{1}{2}\ln\left(\frac{1+x}{1-x} \right)\pm \frac{z_{0,\!975}}{\sqrt{n-3}}=[a,b],$$ où $z_{0,\!975}=1,\!96$. Ensuite en appliquant la fonction inverse qui est $$\mbox{tanh(x)}=\frac{e^{2\,x}-1}{e^{2\,x}+1},$$ on obtient un IC pour $\rho$. Alors, un IC \`a 95\% pour $\rho$ est $$\left[\frac{e^{2\,a}-1}{e^{2\,a}+1},\frac{e^{2\,b}-1}{e^{2\,b}+1} \right]=[0,\!797;0,\!887].$$
  
  Avec R:

  > r<-with(VC,cor(Taille.en.cm,VC))
  > z<-qnorm(0.975)
  > a<-(1/2)*log((1+r)/(1-r))-z/sqrt(n-3)
  > n<-nrow(VC)
  > n<-nrow(VC)
  > r<-with(VC,cor(Taille.en.cm,VC))
  > z<-qnorm(0.975)
  > a<-(1/2)*log((1+r)/(1-r))-z/sqrt(n-3)
  > b<-(1/2)*log((1+r)/(1-r))+z/sqrt(n-3)
  > (exp(2*a)-1)/(exp(2*a)+1)
  [1] 0.7967316
  > (exp(2*b)-1)/(exp(2*b)+1)
  [1] 0.8869721
  

> SAT<-read.csv("sat.csv")
> head(SAT)
         State Expenditure.per.student student.teacher.ratio Avg.teacher.salary Percent.taking verbalSAT99 mathSAT99 totalSAT99
1    "Alabama"                   4.405                  17.2             31.144              8         491       538       1029
2     "Alaska"                   8.963                  17.6             47.951             47         445       489        934
3    "Arizona"                   4.778                  19.3             32.175             27         448       496        944
4   "Arkansas"                   4.459                  17.1             28.934              6         482       523       1005
5 "California"                   4.992                  24.0             41.078             45         417       485        902
6   "Colorado"                   5.443                  18.4             34.571             29         462       518        980
> dim(SAT)
[1] 50  8

> with(SAT,plot(Avg.teacher.salary,totalSAT99,xlab="Salaire Moyen",ylab="SAT (Math et Verbal)"))
> mod<-lm(totalSAT99~Avg.teacher.salary,data=SAT)
> abline(mod)
> summary(mod)

Call:
lm(formula = totalSAT99 ~ Avg.teacher.salary, data = SAT)

Residuals:
     Min       1Q   Median       3Q      Max
-147.125  -45.354    4.073   42.193  125.279

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)
(Intercept)        1158.859     57.659  20.098  < 2e-16 ***
Avg.teacher.salary   -5.540      1.632  -3.394  0.00139 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 67.89 on 48 degrees of freedom
Multiple R-squared:  0.1935,    Adjusted R-squared:  0.1767
F-statistic: 11.52 on 1 and 48 DF,  p-value: 0.001391

> summary(SAT$Percent.taking)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
   4.00    9.00   28.00   35.24   63.00   81.00
> ## definir une variable muette
> SAT$MuetteHaut<-as.numeric(SAT$Percent.taking>=28)

• Si nous avions utilisé la commande SAT$MuetteHaut<-(SAT$Percent.taking>=28), alors nous aurions une variable Booléenne (vrai ou faux) nommé MuetteHaut dans le jeu de données SAT.
• R interprète un vrai comme 1 et un faux comme 0, alors si on force un vecteur Booléen à être numérique, ceci veut dire que TRUE devient 1 et FALSE devient 0.
  
  
  Voici un nuage de points avec la stratification. Qu'observez-vous?
  
  
  
  

  with(SAT,plot(Avg.teacher.salary,totalSAT99,
  xlab="Salaire Moyen",ylab="SAT (Math et Verbal)",
  col=factor(MuetteHaut), cex=1.25))
  with(SAT,text(Avg.teacher.salary,totalSAT99,MuetteHaut,cex=0.75, pos=4))
  

> ## fonction I() nous permet d'appliquer des operations
> ## sur des colonnes du dataframe
> mod<-lm(totalSAT99~MuetteHaut+Avg.teacher.salary+I(MuetteHaut*Avg.teacher.salary),data=SAT)
> mod$coefficients
                       (Intercept)                         MuetteHaut                 Avg.teacher.salary
                      1035.7189251                       -174.6680327                         -0.1733339
I(MuetteHaut * Avg.teacher.salary)
                         1.2390482
> OrdHaut=mod$coefficients[1]+mod$coefficients[2]
> PenteHaut=mod$coefficients[3]+mod$coefficients[4]
> OrdBas=mod$coefficients[1]
> PenteBas=mod$coefficients[3]
> ## parametres (estimes) des deux droites
> c(OrdHaut,PenteHaut,OrdBas,PenteBas)
       (Intercept) Avg.teacher.salary        (Intercept) Avg.teacher.salary
       861.0508923          1.0657143       1035.7189251         -0.1733339

with(SAT,plot(Avg.teacher.salary,totalSAT99,
xlab="Salaire Moyen",ylab="SAT (Math et Verbal)",
col=factor(MuetteHaut), cex=1.25))
with(SAT,text(Avg.teacher.salary,totalSAT99,MuetteHaut,cex=0.75, pos=4))
abline(OrdHaut,PenteHaut)
abline(OrdBas,PenteBas,lty=2)

> mod<-lm(totalSAT99~Avg.teacher.salary+Percent.taking,data=SAT)
> summary(mod)

Call:
lm(formula = totalSAT99 ~ Avg.teacher.salary + Percent.taking,
    data = SAT)

Residuals:
    Min      1Q  Median      3Q     Max
-78.313 -26.731   3.168  18.951  75.590

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)
(Intercept)        987.9005    31.8775  30.991   <2e-16 ***
Avg.teacher.salary   2.1804     1.0291   2.119   0.0394 *
Percent.taking      -2.7787     0.2285 -12.163    4e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 33.69 on 47 degrees of freedom
Multiple R-squared:  0.8056,    Adjusted R-squared:  0.7973
F-statistic: 97.36 on 2 and 47 DF,  p-value: < 2.2e-16

> genou<-read.csv("genou.csv")
> head(genou)
  Temps Groupe  Age
1    29      1 18.3
2    42      1 30.0
3    38      1 26.5
4    40      1 28.1
5    43      1 29.7
6    40      1 27.8
> genou$Groupe<-factor(genou$Groupe)
> levels(genou$Groupe)
[1] "1" "2" "3"

• Le vecteur genou était numérique. Mais en réalité, les nombres 1, 2, et 3 ne sont pas des quantités, mais des étiquettes pour des catégories. Alors, on a utilisé la fonction factor pour convertir le vecteur numérique en facteur (c'est-à-dire une variable catégorique avec R).
• On peut utiliser la fonction constrasts pour afficher le codage pour un facteur. On observe qu'il y a deux variables muettes: une pour identifier le groupe 2 et l'autre pour identifier le groupe 3.

  > contrasts(genou$Groupe)
    2 3
  1 0 0
  2 1 0
  3 0 1
  

• Alors, l'ordonnée à l'origine \(\beta_0=\mu_1\) est la moyenne du group de référence.
• Alors, le coefficient de la variable muette pour le groupe 2 est \(\beta_1=\mu_2-\mu_1\). Ceci est dit l'effet du groupe 2. C'est la différence entre la moyenne du groupe 2 et le groupe de référence.
• Alors, le coefficient de la variable muette pour le groupe 3 est \(\beta_2=\mu_3-\mu_1\). Ceci est dit l'effet du groupe 3. C'est la différence entre la moyenne du groupe 3 et le groupe de référence.
• Si les effets des groupes 2 et 3 sont nuls, alors ceci veut dire que les moyennes sont égaux, et il y a une moyenne en commune pour les trois groupes. Alors, le test de la signification de la régression que \[ H_0:\beta_1=\beta_2=0 ~~ \mbox{ contre } ~~ H_a:\beta_1\neq 0 ~~~ \mbox{ ou } ~~~ \beta_2\neq 0 \] est équivalents aux hypothèses \[ H_0:\mu_1=\mu_2=\mu_3 ~~ \mbox{ contre } ~~ H_a:~ \mbox{les moyennes ne sont pas égaux}. \] Ceci fut la première application d'une ANOVA par R.A. Fisher pour tester l'égalité des moyennes de populations normales et indépendantes.
• Un modèle linéaire qui a seulement des prédicteurs catégoriques est dit un modèle d'ANOVA. Ici nous avons un modèle d'ANOVA à un facteur.

mod<-lm(Temps~Groupe,data=genou)
> summary(mod)

Call:
lm(formula = Temps ~ Groupe, data = genou)

Residuals:
   Min     1Q Median     3Q    Max
  -9.0   -3.0   -0.5    3.0    8.0

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   38.000      1.574  24.149  < 2e-16 ***
Groupe2       -6.000      2.111  -2.842  0.00976 **
Groupe3      -14.000      2.404  -5.824 8.81e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.451 on 21 degrees of freedom
Multiple R-squared:  0.6176,	Adjusted R-squared:  0.5812
F-statistic: 16.96 on 2 and 21 DF,  p-value: 4.129e-05

• On estime que les temps moyens de rétablissement sont \[ \hat{\mu}_1=38 \mbox{ mois}; ~~ \hat{\mu}_2=38-6=32 \mbox{ mois}; ~~ \hat{\mu}_3=38-14=24 \mbox{ mois}. \] L'écart type résiduel est \(s_e=4.451\) mois.
  
  Description: On estime qu'un patient avec un forme inférieure à la moyenne aura en moyenne un temps de rétablissement de 38 mois. Mais, si le patient a une forme moyenne, son temps de rétablissement sera réduite de 6 mois en moyenne. Cette réduction sera 14 mois en moyenne pour les patients ayant un forme supérieure à la moyenne. L'écart type résiduel est de 4,451 mois.
• Avec R, si un modèle linéaire a seulement un prédicteur, alors on peut obtenir le tableau d'ANOVA pour le test de la signifcation de la régression avec la fonction anova. Ceci est le test de l'égalité des moyennes. On conclue qu'il y a des différences entre les temps moyen de rétablissement selon la forme du patient \((F(2,21)=16.96; p<0.0001)\).

  > anova(mod)
  Analysis of Variance Table
  
  Response: Temps
            Df Sum Sq Mean Sq F value    Pr(>F)
  Groupe     2    672  336.00  16.962 4.129e-05 ***
  Residuals 21    416   19.81
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  

• On peut aussi obtenir le test de la signification de la régression en utilisant un test linéaire générale. On compare le modèle complet avec le modèle réduit est le modèle \(E\{Y\}=\beta_0\).

  > # modèle réduit
  > mod0<-lm(Temps~1,data=genou)
  > anova(mod0,mod)
  Analysis of Variance Table
  
  Model 1: Temps ~ 1
  Model 2: Temps ~ Groupe
    Res.Df  RSS Df Sum of Sq      F    Pr(>F)
  1     23 1088
  2     21  416  2       672 16.962 4.129e-05 ***
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  

• Ici nous avons analysé des données d'une étude observationelle. Il est possible que la différence entre les groupes que nous avons observé peut être expliqué par une variable confusionelle. C'est une pratique commune des études bio-médicales d'ajuster pour l'âge du patient et le sexe du patient. Dans la plupart des cas, l'effet du sexe est très important. Il est aussi une pratique commune de séparer les sexes dans les applications médicales. C'est-à-dire, une étude va souvent seulement avoir des hommes ou seulement des femmes. Les chercheurs de cette étude veulent ajuster pour l'âge. Il est possible que l'âge des patients n'est pas équitablement distribué dans les groupes de la forme. ll est possible qu'on ait plus de patients âgés dans le groupe de forme inférieure à la moyenne, et que c'est ceci qui est la raison qu'on observe des différences entre les groupes. Pour ajuster pour l'âge, on ajoute l'âge comme un prédicteur quantitatif au modèle.

  > mod<-lm(Temps~Groupe+Age,data=genou)
  > summary(mod)
  
  Call:
  lm(formula = Temps ~ Groupe + Age, data = genou)
  
  Residuals:
       Min       1Q   Median       3Q      Max
  -1.03891 -0.36892  0.05891  0.33098  0.89991
  
  Coefficients:
              Estimate Std. Error t value Pr(>|t|)
  (Intercept)  7.43169    0.86106   8.631 3.54e-08 ***
  Groupe2     -1.84738    0.28694  -6.438 2.80e-06 ***
  Groupe3     -8.72289    0.33296 -26.198  < 2e-16 ***
  Age          1.16729    0.03201  36.461  < 2e-16 ***
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  Residual standard error: 0.5552 on 20 degrees of freedom
  Multiple R-squared:  0.9943,	Adjusted R-squared:  0.9935
  F-statistic:  1170 on 3 and 20 DF,  p-value: < 2.2e-16
   

  
  
  Remarques:
  • Un modèle linéaire avec des prédicteurs quantitatifs et catégoriques est dit un modèle linéaire général. Ici nous avons un modèle linéaire général d'ordre 1 qui décrit le temps de rétablissement selon la forme pré-chirurgie du patient et son âge.
  • Le modèle est significatif \(F(3,20)=1170; p<0,\!0001\); \(R^2=0,\!9943\) et l'erreur type résiduelle est 0,5552 mois.
  • Le modèle estimée est \[ \begin{eqnarray*} \hat{\mbox{temps}} &=& 7,\!43169-1,\!84738\,\,I\{\mbox{Groupe}=2\}-8,\!72289\,I\{\mbox{Groupe}=3\}+1,\!16729\,\mbox{âge}. \\ \end{eqnarray*} \]
  • L'estimation de l'ordonnée à l'origine est $b_0=7,\!43169$. Ceci est le temps de rétablissement pour un patient dans le groupe 1, mais avec un âge de 0. Ceci n'est pas signifiant. Pour donner un sens à l'ordonnée à l'origine, on peut centrer l'âge autour de sa moyenne. L'âge moyen des participants est 23,757 ans. On pourrait centrer autour d'un âge de 24 pour simplifier la description.

    > mean(genou$Age)
    [1] 23.575
    

  • On ajoute au jeu de données une variable qui sera l'âge centré à 24.

      > genou$Age.c<-genou$Age-24
     

  • On ajuste le modèle de nouveau, mais plutôt avec la variable âge centrée.

     > mod<-lm(Temps~Groupe+Age.c,data=genou)
    > summary(mod)
    
    Call:
    lm(formula = Temps ~ Groupe + Age.c, data = genou)
    
    Residuals:
         Min       1Q   Median       3Q      Max
    -1.03891 -0.36892  0.05891  0.33098  0.89991
    
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)
    (Intercept) 35.44656    0.20842 170.070  < 2e-16 ***
    Groupe2     -1.84738    0.28694  -6.438  2.8e-06 ***
    Groupe3     -8.72289    0.33296 -26.198  < 2e-16 ***
    Age.c        1.16729    0.03201  36.461  < 2e-16 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 0.5552 on 20 degrees of freedom
    Multiple R-squared:  0.9943,	Adjusted R-squared:  0.9935
    F-statistic:  1170 on 3 and 20 DF,  p-value: < 2.2e-16
     

    Le modèle estimé est \[ \begin{eqnarray*} \hat{\mbox{temps}} &=& 35,\!44656-1,\!84738\,\,I\{\mbox{Groupe}=2\}-8,\!72289\,I\{\mbox{Groupe}=3\}+1,\!16729\,(\mbox{âge}-24). \\ &=&\left\{ \begin{array}{ll} 35,\!44656+1,\!16729\,(\mbox{âge}-24), & \mbox{Groupe}=1\\ 33,\!59918+1,\!16729\,(\mbox{âge}-24) & \mbox{Groupe}=2\\ 26,\,72367+1,\!16729\,(\mbox{âge}-24) & \mbox{Groupe}=3\\ \end{array} \right. \end{eqnarray*} \]
    
    Pour un patient âgé de 24 ans de forme inférieure à la moyenne, on estime que le temps moyen de rétablissement est 35,4 mois. Si ce patient de 24 ans a une forme moyenne, ceci réduit le temps moyen de rétablissement 1,8 mois, et cette réduction est de 8,7 mois pour un patient de 24 ans avec une forme supérieure à la moyenne. Pour chaque année qu'on ajoute à l'âge du patient, on estime que le temps moyen de rétablissement va augmenter de 1,2 mois. L'écart type résiduelle est 0,552 mois.
  
  
  
  

  Exemple 24 : Test de la signification de la régression

  
  
  
  Considérons un modèle linéaire avec la fonction de la moyenne suivante $$E\{Y\}=\beta_0+\beta_1\,x_1+\ldots+\beta_{p-1}\,x_{p-1}.$$
  
  Une façon de décrire la valeur explicative du modèle est de vérifier que le modèle est significatif, dans le sense qu'il y a au moins une variable explicative significative.
  
  Le test de la signification de la régression est de confronter $$H_0~:~ \beta_j=0, ~~ \mbox{ pour } ~ j=1,2,\ldots,p-1 ~~~~ \mbox{ contre } ~~~~ H_a~:~ \beta_j\neq 0, ~~ \mbox{ pour au moins un $j$ }.$$
  
  La statistique du test est $$F^*=\frac{MSR}{MSE}=\frac{SSR/(p-1)}{SSE/(n-p)},$$ qui suit une loi $F(p-1,n-p)$ si $H_0$ est vraie. Des grandes valeurs de $F^*$ sont considéré comme des preuves contre $H_0$.
  
  Exemple [Dwayne Studios]: Nous allons considérer les données dans le fichier Dwayne.csv. Dwayne studios est une compagnie qui a des studios dans $n=21$ villes de taille moyenne. Leur spécialité est la photographie des enfants. Les variables sont des descriptions des studios dans chacune de ces villes. Elles sont les ventes $y$, le nombre d'enfants d'\^ages 16 ou moins dans la communaut\'e $x_1$, et le revenu disponible par habitant $x_2$ (en milliers de dollars).
  
  Le mod\`ele de r\'egression estim\'ee est $$\hat{y}=-68,\!8571+0,\!00145\,x_1+9,\!3655\,x_2.$$

  > dwayne<-read.csv("Dwayne.csv")
  > head(dwayne)
    sales num.children disposible.income
  1 174.4        68500              16.7
  2 164.4        45200              16.8
  3 244.2        91300              18.2
  4 154.6        47800              16.3
  5 181.6        46900              17.3
  6 207.5        66100              18.2
  > mod<-lm(sales~num.children+disposible.income,data=dwayne)
  > coefficients(mod)
        (Intercept)      num.children disposible.income
       -68.85707315        0.00145456        9.36550038
  

  Est-ce qu'on peut utiliser $x_1$ et $x_2$ pour pr\'evoir les ventes dans une autre communit\'e? En d'autres-mots, est-ce que la régression est significative?
  
  On veut tester $H_0:\beta_1=\beta_2=0$ contre $H_a:$ au moins un des coefficients $\beta_1, \beta_2$ n'est pas zéro.
  
  La statistique du test est $$F^*=\frac{SSR/(p-1)}{SSE/(n-p)}=99,\!1$$ et la valeur-$p$ est $P(F(2,18)>99,\!1)=1,\!921\times 10^{-10}.$ La régression est significative.

  > summary(mod)
  
  Call:
  lm(formula = sales ~ num.children + disposible.income, data = dwayne)
  
  Residuals:
       Min       1Q   Median       3Q      Max
  -18.4239  -6.2161   0.7449   9.4356  20.2151
  
  Coefficients:
                      Estimate Std. Error t value Pr(>|t|)
  (Intercept)       -6.886e+01  6.002e+01  -1.147   0.2663
  num.children       1.455e-03  2.118e-04   6.868    2e-06 ***
  disposible.income  9.366e+00  4.064e+00   2.305   0.0333 *
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  Residual standard error: 11.01 on 18 degrees of freedom
  Multiple R-squared:  0.9167,    Adjusted R-squared:  0.9075
  F-statistic:  99.1 on 2 and 18 DF,  p-value: 1.921e-10
  

  
  
  Tester pour la signication de la régression avec le test linéaire général: On veut tester $H_0:\beta_1=\beta_2=0$. Ceci est équivalent à confronter $$H_0: E\{Y\}=\beta_0 ~~~ \mbox{ contre } ~~~ H_a: E\{Y\}=\beta_0+\beta_1\,x_1+\beta_2\,x_2.$$ Nous ajustons les deux modèles et on les compare avec la fonction anova. La statistique du test est $$F^*=\frac{\mbox{Extra}/(p-q)}{\mbox{SSE}/(n-p)}=99,\!103,$$ et la valeur-$p$ est $P(F(2,18)>99,\!103)=1,\!92\times 10^{-10}$.

  > mod<-lm(sales~num.children+disposible.income,data=dwayne)
  > mod0<-lm(sales~1,data=dwayne)
  > anova(mod0,mod)
  Analysis of Variance Table
  
  Model 1: sales ~ 1
  Model 2: sales ~ num.children + disposible.income
    Res.Df     RSS Df Sum of Sq      F    Pr(>F)
  1     20 26196.2
  2     18  2180.9  2     24015 99.103 1.921e-10 ***
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  

  
  
  
  Commentaires: Les deux approches pour le test de la signification de la régression sont équivalentes.
  
  
  
  
  
  
  

  Exemple 25 : La loi F non-centrée

  
  
  
  Considérons un test linéaire général avec $p=7$, $q=2$ et $n=20$. La statistique du test $F^*$ suit une loi $F(p-q=5, n-p=13, \lambda)$, o\`u $$\lambda=\frac{E\{\mathbf{Y}\}'(H-H_0)E\{\mathbf{Y}\}}{\sigma^2}=n\,f^2,$$ et $H$ est la matrice chapeau du mod\`ele complet et $H_0$ est la matrice chapeau du mod\`ele r\'eduit, et $f^2$ est le $f^2$ de Cohen.
  
  Voici un graphe de la densité de probabilité de la loi $F(p-q=5, n-p=13, \lambda)$, pour $\lambda=0$, $\lambda=1,$\lambda=2$, et$\lambda=3$. Nous avons aussi superimposé une droite verticale à la valeur critique$F(0,\!95;5,13)$. Rappel: On rejette$H_0$si$F^*>F(0,\!95;5,13)$. Plus que$\lambda=n\,f^2$est grand, plus que chance de rejetter$H_0$ augmente.
  
  
  
  
  
  

  curve(df(x,5,13,0),0,6,lty=1,ylab="densité",xlab="F")
  curve(df(x,5,13,1),0,6,lty=2,add=TRUE,col="red")
  curve(df(x,5,13,2),0,6,lty=3,add=TRUE,col="blue")
  curve(df(x,5,13,3),0,6,lty=4,add=TRUE,col="purple")
  abline(v=qf(0.95,5,13))
  legend(4, 0.5, legend=c("n.c.=0","n.c.=1","n.c.=2","n.c.=3"),
         col=c("black","red", "blue","purple"), lty=1:4,
         title="Non Centralité")
  

  
  
  Calcul de puissance: La puissance d'un test est définie comme étant la probabilité de rejetter $H_0$. Quand $H_0$ est vraie, la puissance est $P(\mbox{rejet $H_0$}|\mbox{$H_0$ est vraie})=\alpha$. Mais, si $H_a$ est vraie la puissance sera un fonction de la taille de l'échantillon et de la taille de l'effet. Dans le cas du test linéaire générale, la puissance est une fonction du paramètre de non-centralité $\lambda=n\,f^2$, où $f^2$ est le $f^2$ de Cohen. La puissance du test linéaire générale est $$\mbox{puissance}=P(F(p-q,n-p,\lambda=n\,f^2)>F(0,95;p-q,n-p)).$$
  
  Ecrivons notre propre fonction qui calcul la puissance avec R.

  power.partiel.F <- function(p,q,n,f2,alpha)
  {
     1-pf(qf(1-alpha,p-q,n-p),p-q,n-p,n*f2)
  }
  

  
  
  Exemple: Considérons un test linéaire générale avec $p=7$, $q=2$ et $n=20$. Calculons la puissance pour (a) $f^2=0,\!1$; (b) $f^2=2,\!5$.
  
  (a) La puissance est $$\mbox{puissance}=P(F(p-q=5,n-p=13,\lambda=n\,f^2=2)>F(0,95;p-q=5,n-p=13))=P(F(5,13,2)>3.025438)=0.1194.$$

  > qf(0.95,5,13)
  [1] 3.025438
  > power.partiel.F(p=7,q=2,n=20,f2=0.1,alpha=0.05)
  [1] 0.119406
  

  
  
  (b) La puissance est $$\mbox{puissance}=P(F(p-q=5,n-p=13,\lambda=n\,f^2=50)>F(0,95;p-q=5,n-p=13))=P(F(5,13,50)>3.025438)=0.9975.$$ ]

  >  power.partiel.F(p=7,q=2,n=20,f2=2.5,alpha=0.05)
  [1] 0.9975355
  

  
  
  Calcul d'une taille d'échantillon: (a) Nous voulons planifier une étude. Nous allons utiliser un test linéaire générale avec $p=7$, $q=2$. Nous voulons déterminer la taille de l'échantillon $n$ adéquate afin d'identifier une taille de l'effet de $f^2=1,\!5$ avec une puissance de 80\% \`a un niveau de signification de $\alpha=5\%$.
  
  Puisque $n> p$, alors nous allons évaluer la puissance du test à $n=8, 8, 9, 10,$ et ainsi de suite jusqu\'a ce que nous avons une puissance de 80%. On observe que la taille de l'échantillon requise est $n=16$ et la puissance correspondante est $81,\!88\%$.
  
  

  > p<-7; q<-2; alpha<-0.05; f2<-1.5
  > n<-p
  > power<-0
  > while (power <0.8)
  + {
  + n<-n+1
  + power<-power.partiel.F(p,q,n,f2,alpha)
  + }
  > power
  [1] 0.8187186
  > n
  [1] 16
  

  
  
  (b) Vérifions qu'à $n=15$, la puisance sera plus petite que 80%. La puissance à $n=15$, si $f^2=2,\!5$, est $$\mbox{puissance}=P(F(p-q=5,n-p=15-7=8,\lambda=n\,f^2=(15)(1,\!5)=22,\!5)>F(0,\!95;p-q=5,n-p=8))=P(F(5,8,22,\!5)> 3,\!687499)=0.7594.$$

  > qf(0.95,5,8)
  [1] 3.687499
  > power.partiel.F(p=7,q=2,n=15,f2=1.5,alpha=0.05)
  [1] 0.7593985
  

  
  
  
  
  
  
  
  
  
  
  
  
  
  

  Exemple 23 : Comparer des modèles emboités pour tester des hypothèses (Test F partiel)

  
  
  
  Considérons un modèle linéaire avec la matrice du plan $X$ de taille $n\times p$. Supposons que nous voulons tester $H_0$ qui nous donne un modèle réduit qui est emboité dans notre modèle. Supposons que $X_0$ de taille $n\times q$ est la matrice du plan du modèle réduit, où les colonnes de $X_0$ appartiennent à l'espace colonne de $X$. On peut tester $H_0$ avec la statistique $F$ suivante: $$F^*=\frac{\mbox{ExtraSS}/(p-q)}{\mbox{SSE}/(n-p)}, ~~ \mbox{ où } ~~ \mbox{ExtraSS}=\mbox{SSE}(R)-\mbox{SSE}.$$ Plus que ExtraSS est grande, plus que nous avons des preuves contre $H_0$. Sous $H_0$, la statistique suit une loi $F(p-q,n-p)$. Alors, la valeur-$p$ est $P(F(p-q,n-p)>F^*)$, où $F^*$ est la valeur observ\'ee de la statistique.
  
  Exemple [Rétablissement]: Considérons les données d'une étude observationnelle bio-médicale. Les participants de cette étude dans une étude a subi une opération du genou. La variable dépendante est le temps de rééducation, et la variable explicative est une variable catégorique décrivant la forme du patient avant une intervention chirurgicale (1 = inférieur à la moyenne; 2 = moyen; 3 = supérieur à la moyenne). Les 24 participants sont des hommes de 18 ans à 30 ans. Nous avons aussi l'âge du participant dans le jeu de données. Les données sont dans le fichier : genou.csv. On importe les données et on affiche quelques rangés du jeu de données. On change type du vecteur \verb|genou| en vecteur catégorique et on affiche ces niveaux.
  
  

  > genou<-read.csv("genou.csv")
  > head(genou)
    Temps Groupe  Age
  1    29      1 18.3
  2    42      1 30.0
  3    38      1 26.5
  4    40      1 28.1
  5    43      1 29.7
  6    40      1 27.8
  > genou$Groupe<-factor(genou$Groupe)
  > levels(genou$Groupe)
  [1] "1" "2" "3"
  

  
  
  On veut décrire le temps de rétablissement selon la forme, mais en ajustant l'âge. On commence par tester qu'il n'y pas d'interactions entre l'âge et la forme.
  
  Nous allons ajuster le modèle linéaire suivant: $$E\{Y\}=\beta_0+\beta_1\,\,I\{\mbox{Groupe}=2\} +\beta_2\,I\{\mbox{Groupe}=3\}+ \beta_3\,(\mbox{Age}) +\beta_4\,(\mbox{Age})\,I\{\mbox{Groupe}=2\}+\beta_5\,(\mbox{Age})\,I\{\mbox{Groupe}=3\}.$$
  
  Nous voulons tester l'hypothèse nulle qu'il n'y a pas d'interactions entre la forme et l'âge. $$H_0~:~ E\{Y\}=\beta_0+\beta_1\,\,I\{\mbox{Groupe}=2\} +\beta_2\,I\{\mbox{Groupe}=3\}+ \beta_3\,(\mbox{Age})$$ contre $$H_a~:~E\{Y\}=\beta_0+\beta_1\,\,I\{\mbox{Groupe}=2\} +\beta_2\,I\{\mbox{Groupe}=3\}+ \beta_3\,(\mbox{Age}) +\beta_4\,(\mbox{Age})\,I\{\mbox{Groupe}=2\}+\beta_5\,(\mbox{Age})\,I\{\mbox{Groupe}=3\}.$$
  
  On ajuste les deux modèles et on les compare avec la fonction anova. N.B. En imposant les contraintes de $H_0$, on a perdu deux paramètres. Alors, le nombre de degrés de liberté du numérateur est 2.
  
  

  > mod<-lm(formula = Temps ~ Groupe + Age+Groupe*Age, data = genou)
  > mod<-lm(formula = Temps ~ Groupe + Age+Groupe*Age, data = genou)
  > mod0<-lm(formula = Temps ~ Groupe + Age, data = genou)
  > anova(mod0,mod)
  Analysis of Variance Table
  
  Model 1: Temps ~ Groupe + Age
  Model 2: Temps ~ Groupe + Age + Groupe * Age
    Res.Df    RSS Df Sum of Sq      F Pr(>F)
  1     20 6.1657
  2     18 5.9439  2   0.22184 0.3359 0.7191
  

  La statistique du test est $$F^*=\frac{\mbox{ExtraSS}/(p-q)}{SSE/(n-p)}=\frac{0,\!22184/2}{5,\!9439/18}=0.3359.$$ La valeur-$p$=P(F(2,18)>0.3359)=0.7191)$. Les preuves qu'il y ait des interactions entre l'âge et la forme ne sont pas significatives.
  
  
  
  Signification des variables explicatives d'un modèle additif: Nous allons utiliser un modèle linéaire additif pour décrire le lien entre le temps de rétablissements et les variables explicatives qui sont la forme pré-opération et l'âge.
  
  La fonction de la moyenne du modèle est $$E\{Y\}=\beta_0+\beta_1\,\,I\{\mbox{Groupe}=2\} +\beta_2\,I\{\mbox{Groupe}=3\}+ \beta_3\,(\mbox{Age}).$$
  
  
  • Test pour la signification de l'âge: On veut tester $H_0:\beta_3=0$ contre $H_a:\beta_3\neq 0$. On peut utiliser un statistique $t$ ou $F$ pour tester ces hypothèses.
    
    Avec une statistique $t$~: La valeur observée de la statistique est $t^*=b_3/s\{b3\}=36,\!4608.$ La valeur-$p$ est $2\,P(T(20)>36,\!4608)=9,1\times 10^{-20}$. L'âge est une variable explicative significative du temps de rétablissement.

    > mod<-lm(formula = Temps ~ Groupe + Age, data = genou)
    > summary(mod)$coeff
                 Estimate Std. Error    t value     Pr(>|t|)
    (Intercept)  7.431688 0.86106382   8.630821 3.538937e-08
    Groupe2     -1.847379 0.28694289  -6.438141 2.802274e-06
    Groupe3     -8.722893 0.33296397 -26.197708 5.914908e-17
    Age          1.167286 0.03201483  36.460804 9.075467e-20
    > mod$df.residual
    [1] 20
    

    
    
    Avec une statistique $F$~: On veut tester $$H_0~:~E\{Y\}=\beta_0+\beta_1\,\,I\{\mbox{Groupe}=2\} +\beta_2\,I\{\mbox{Groupe}=3\} ~~~ \mbox{ contre } ~~~ H_a~:~E\{Y\}=\beta_0+\beta_1\,\,I\{\mbox{Groupe}=2\} +\beta_2\,I\{\mbox{Groupe}=3\}+ \beta_3\,(\mbox{Age}).$$
    
    La statistique du test est $$F^*=\frac{\mbox{ExtraSS}/(p-q)}{\mbox{SSE}/(n-p)}=1329,\!4.$$ La valeur-$p$ du test est $P(F(1,20)>1329,\!4)=2,\!2\times 10^{-16}$.
    
    

    > mod<-lm(formula = Temps ~ Groupe + Age, data = genou)
    > mod0<-lm(formula = Temps ~ Groupe , data = genou)
    > anova(mod0,mod)
    Analysis of Variance Table
    
    Model 1: Temps ~ Groupe
    Model 2: Temps ~ Groupe + Age
      Res.Df    RSS Df Sum of Sq      F    Pr(>F)
    1     21 416.00
    2     20   6.17  1    409.83 1329.4 < 2.2e-16 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    

    
    
    Remarque: Les deux tests sont équivalents. On peut démontrer que $(t^*)^2=F^*$.
    
    
  • Test pour la signification de la forme pré-opération~: La fonction de la moyenne du modèle est $$E\{Y\}=\beta_0+\beta_1\,\,I\{\mbox{Groupe}=2\} +\beta_2\,I\{\mbox{Groupe}=3\}+ \beta_3\,(\mbox{Age}).$$ On veut tester $$H_0:\beta_1=0, \beta_2=0 ~~~ \mbox{contre} ~~~ H_a: \mbox{$H_0$ est fausse}.$$ Ici on ne peut pas utiliser un test $t$. Mais, on peut utiliser un test $F$. On veut tester $$H_0~:~E\{Y\}=\beta_0+\beta_3\,(\mbox{Age}) ~~~ \mbox{ contre } ~~~ H_a~:~E\{Y\}=\beta_0+\beta_1\,\,I\{\mbox{Groupe}=2\} +\beta_2\,I\{\mbox{Groupe}=3\}+ \beta_3\,(\mbox{Age}).$$ La statistique du test est $$F^*=\frac{\mbox{ExtraSS}/(p-q)}{\mbox{SSE}/(n-p)}=399,\!11.$$ La valeur-$p$ du test est $P(F(2,20)>399,\!11)=2,\!2\times 10^{-16}$. La forme pré-opération est une variable explicative significative.

    > mod<-lm(formula = Temps ~ Groupe + Age, data = genou)
    > mod0<-lm(formula = Temps ~ Age, data = genou)
    > anova(mod0,mod)
    Analysis of Variance Table
    
    Model 1: Temps ~ Age
    Model 2: Temps ~ Groupe + Age
      Res.Df     RSS Df Sum of Sq      F    Pr(>F)
    1     22 252.249
    2     20   6.166  2    246.08 399.11 < 2.2e-16 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    

  
  
  
  
  Obtenir la signification de toutes les variables explicatives d'un modèle: On peut utiliser la fonction anova, mais c'est pas exactement ce que nous allons vouloir. Ensuite, nous allons vous montrer la fonction drop1 et la fonction Anova du package car.
  • La fonction anova sur un seul objet lm nous donne des tests séquentielles. Remarquer que si Age est la deuxième variable explicative, alors son $F^*$ est 1329. Tandis que si Age est la première variable explicative, son $F^*$ est 399. Lequel est correct?

    > mod<-lm(formula = Temps ~ Groupe + Age, data = genou)
    > anova(mod)
    Analysis of Variance Table
    
    Response: Temps
              Df Sum Sq Mean Sq F value    Pr(>F)
    Groupe     2 672.00  336.00  1089.9 < 2.2e-16 ***
    Age        1 409.83  409.83  1329.4 < 2.2e-16 ***
    Residuals 20   6.17    0.31
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    > mod<-lm(formula = Temps ~ Age+Groupe, data = genou)
    > anova(mod)
    Analysis of Variance Table
    
    Response: Temps
              Df Sum Sq Mean Sq F value    Pr(>F)
    Age        1 835.75  835.75 2710.95 < 2.2e-16 ***
    Groupe     2 246.08  123.04  399.11 < 2.2e-16 ***
    Residuals 20   6.17    0.31
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    

    Si on utilise Temps ~ Groupe + Age, alors les hypothèses pour Groupe sont $$H_0~:~E\{Y\}=\beta_0 ~~~ \mbox{ contre } ~~~ H_a~:~E\{Y\}=\beta_0+\beta_1\,\,I\{\mbox{Groupe}=2\} +\beta_2\,I\{\mbox{Groupe}=3\}.$$ et les hypothèses pour Age sont $$H_0~:~E\{Y\}=\beta_0+\beta_1\,\,I\{\mbox{Groupe}=2\} +\beta_2\,I\{\mbox{Groupe}=3\}. ~~~ \mbox{ contre } ~~~ H_a~:~E\{Y\}=\beta_0+\beta_1\,\,I\{\mbox{Groupe}=2\} +\beta_2\,I\{\mbox{Groupe}=3\}+ \beta_3\,(\mbox{Age}).$$ Alors, seulement le test pour Age est correct, dans le sense qu'on ajuste pour Groupe.
    
    Si on utilise Temps ~ Age + Groupe, alors les hypothèses pour Age sont $$H_0~:~E\{Y\}=\beta_0. ~~~ \mbox{ contre } ~~~ H_a~:~E\{Y\}=\beta_0+ \beta_1\,(\mbox{Age}).$$ et les hypothèses pour Groupe sont $$H_0~:~E\{Y\}=\beta_0+\beta_1\,(\mbox{Age}). ~~~ \mbox{ contre } ~~~ H_a~:~E\{Y\}=\beta_0+\beta_1\,(\mbox{Age})+ \beta_2\,\,I\{\mbox{Groupe}=2\} +\beta_3\,I\{\mbox{Groupe}=3\}.$$ Alors, seulement le test pour Groupe est correct, dans le sense qu'on ajuste pour Age.
    
    
  • Avec la fonction drop1, on peut tester pour la signification de chacune des variables explicatives en ajustant pour les autres variables dans le modèle.

    > mod<-lm(formula = Temps ~ Groupe + Age, data = genou)
    > drop1(mod, test="F")
    Single term deletions
    
    Model:
    Temps ~ Groupe + Age
           Df Sum of Sq    RSS     AIC F value    Pr(>F)
                    6.17 -24.617
    Groupe  2    246.08 252.25  60.457  399.11 < 2.2e-16 ***
    Age     1    409.83 416.00  74.463 1329.39 < 2.2e-16 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    

  
  
  
  

  Exemple 24 : Test de la signification de la régression

  
  
  
  Considérons un modèle linéaire avec la fonction de la moyenne suivante $$E\{Y\}=\beta_0+\beta_1\,x_1+\ldots+\beta_{p-1}\,x_{p-1}.$$
  
  Une façon de décrire la valeur explicative du modèle est de vérifier que le modèle est significatif, dans le sense qu'il y a au moins une variable explicative significative.
  
  Le test de la signification de la régression est de confronter $$H_0~:~ \beta_j=0, ~~ \mbox{ pour } ~ j=1,2,\ldots,p-1 ~~~~ \mbox{ contre } ~~~~ H_a~:~ \beta_j\neq 0, ~~ \mbox{ pour au moins un $j$ }.$$
  
  La statistique du test est $$F^*=\frac{MSR}{MSE}=\frac{SSR/(p-1)}{SSE/(n-p)},$$ qui suit une loi $F(p-1,n-p)$ si $H_0$ est vraie. Des grandes valeurs de $F^*$ sont considéré comme des preuves contre $H_0$.
  
  Exemple [Dwayne Studios]: Nous allons considérer les données dans le fichier Dwayne.csv. Dwayne studios est une compagnie qui a des studios dans $n=21$ villes de taille moyenne. Leur spécialité est la photographie des enfants. Les variables sont des descriptions des studios dans chacune de ces villes. Elles sont les ventes $y$, le nombre d'enfants d'\^ages 16 ou moins dans la communaut\'e $x_1$, et le revenu disponible par habitant $x_2$ (en milliers de dollars).
  
  Le mod\`ele de r\'egression estim\'ee est $$\hat{y}=-68,\!8571+0,\!00145\,x_1+9,\!3655\,x_2.$$

  > dwayne<-read.csv("Dwayne.csv")
  > head(dwayne)
    sales num.children disposible.income
  1 174.4        68500              16.7
  2 164.4        45200              16.8
  3 244.2        91300              18.2
  4 154.6        47800              16.3
  5 181.6        46900              17.3
  6 207.5        66100              18.2
  > mod<-lm(sales~num.children+disposible.income,data=dwayne)
  > coefficients(mod)
        (Intercept)      num.children disposible.income
       -68.85707315        0.00145456        9.36550038
  

  Est-ce qu'on peut utiliser $x_1$ et $x_2$ pour pr\'evoir les ventes dans une autre communit\'e? En d'autres-mots, est-ce que la régression est significative?
  
  On veut tester $H_0:\beta_1=\beta_2=0$ contre $H_a:$ au moins un des coefficients $\beta_1, \beta_2$ n'est pas zéro.
  
  La statistique du test est $$F^*=\frac{SSR/(p-1)}{SSE/(n-p)}=99,\!1$$ et la valeur-$p$ est $P(F(2,18)>99,\!1)=1,\!921\times 10^{-10}.$ La régression est significative.

  > summary(mod)
  
  Call:
  lm(formula = sales ~ num.children + disposible.income, data = dwayne)
  
  Residuals:
       Min       1Q   Median       3Q      Max
  -18.4239  -6.2161   0.7449   9.4356  20.2151
  
  Coefficients:
                      Estimate Std. Error t value Pr(>|t|)
  (Intercept)       -6.886e+01  6.002e+01  -1.147   0.2663
  num.children       1.455e-03  2.118e-04   6.868    2e-06 ***
  disposible.income  9.366e+00  4.064e+00   2.305   0.0333 *
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  Residual standard error: 11.01 on 18 degrees of freedom
  Multiple R-squared:  0.9167,    Adjusted R-squared:  0.9075
  F-statistic:  99.1 on 2 and 18 DF,  p-value: 1.921e-10
  

  
  
  Tester pour la signication de la régression avec le test linéaire général: On veut tester $H_0:\beta_1=\beta_2=0$. Ceci est équivalent à confronter $$H_0: E\{Y\}=\beta_0 ~~~ \mbox{ contre } ~~~ H_a: E\{Y\}=\beta_0+\beta_1\,x_1+\beta_2\,x_2.$$ Nous ajustons les deux modèles et on les compare avec la fonction anova. La statistique du test est $$F^*=\frac{\mbox{Extra}/(p-q)}{\mbox{SSE}/(n-p)}=99,\!103,$$ et la valeur-$p$ est $P(F(2,18)>99,\!103)=1,\!92\times 10^{-10}$.

  > mod<-lm(sales~num.children+disposible.income,data=dwayne)
  > mod0<-lm(sales~1,data=dwayne)
  > anova(mod0,mod)
  Analysis of Variance Table
  
  Model 1: sales ~ 1
  Model 2: sales ~ num.children + disposible.income
    Res.Df     RSS Df Sum of Sq      F    Pr(>F)
  1     20 26196.2
  2     18  2180.9  2     24015 99.103 1.921e-10 ***
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  

  
  
  
  Commentaires: Les deux approches pour le test de la signification de la régression sont équivalentes.
  
  
  
  
  
  
  

  Exemple 25 : La loi F non-centrée

  
  
  
  Considérons un test linéaire général avec $p=7$, $q=2$ et $n=20$. La statistique du test $F^*$ suit une loi $F(p-q=5, n-p=13, \lambda)$, o\`u $$\lambda=\frac{E\{\mathbf{Y}\}'(H-H_0)E\{\mathbf{Y}\}}{\sigma^2}=n\,f^2,$$ et $H$ est la matrice chapeau du mod\`ele complet et $H_0$ est la matrice chapeau du mod\`ele r\'eduit, et $f^2$ est le $f^2$ de Cohen.
  
  Voici un graphe de la densité de probabilité de la loi $F(p-q=5, n-p=13, \lambda)$, pour $\lambda=0$, $\lambda=1,$\lambda=2$, et$\lambda=3$. Nous avons aussi superimposé une droite verticale à la valeur critique$F(0,\!95;5,13)$. Rappel: On rejette$H_0$si$F^*>F(0,\!95;5,13)$. Plus que$\lambda=n\,f^2$est grand, plus que chance de rejetter$H_0$ augmente.
  
  
  
  
  
  

  curve(df(x,5,13,0),0,6,lty=1,ylab="densité",xlab="F")
  curve(df(x,5,13,1),0,6,lty=2,add=TRUE,col="red")
  curve(df(x,5,13,2),0,6,lty=3,add=TRUE,col="blue")
  curve(df(x,5,13,3),0,6,lty=4,add=TRUE,col="purple")
  abline(v=qf(0.95,5,13))
  legend(4, 0.5, legend=c("n.c.=0","n.c.=1","n.c.=2","n.c.=3"),
         col=c("black","red", "blue","purple"), lty=1:4,
         title="Non Centralité")
  

  
  
  Calcul de puissance: La puissance d'un test est définie comme étant la probabilité de rejetter $H_0$. Quand $H_0$ est vraie, la puissance est $P(\mbox{rejet $H_0$}|\mbox{$H_0$ est vraie})=\alpha$. Mais, si $H_a$ est vraie la puissance sera un fonction de la taille de l'échantillon et de la taille de l'effet. Dans le cas du test linéaire générale, la puissance est une fonction du paramètre de non-centralité λ=n