
MAT 3775 - Analyse de la régression
Table des matières
Télécharger et installer R
Exemple 1 : Travailler à la console
Exemple 2 : Affection d'une variable
Exemple 3 : L'éditeur de R.
Exemple 4 : Statistiques descriptives (univariées)
Exemple 5 : Statistiques descriptives (bivariées)
Exemple 6 : Loi de probabilité
Exemple 7 : Importer des données d'un fichier
Exemple 8 : Nuage de points et la surperposition d'une courbe lisse
Exemple 9 : Opérations sur une colonne
Exemple 10 : Ajuster un modèle linéaire avec la fonction lm
Exemple 12 : L'analyse de variance avec un objet lm
Exemple 13 : Un test F partiel pour comparer deux modèles
Exemple 14 : Une variable explicative binaire
Exemple 15 : La signification d'une corrélation
Exemple 21 : Introduction au modèle de régression linéaire multiple (aussi dit modèle linéaire)
Exemple 23 : Comparer des modèles emboités pour tester des hypothèses (Test F partiel)
Exemple 24 : Test de la signifcation de la régression
Exemple 25 : La loi F non-centrée et la puissance
Exemple 31 : La régression polynomiale
Exemple 32 : Une spline linéaire
Exemple 33 : Interaction entre prédicteurs
Exemple 34 : Test de l'inadéquation de l'ajustement (Test for lack-of-fit)
Exemple 35 : Diagnostique
Exemple 36 : Diagnostique pour l'homoscédascité
Exemple 37 : Correction de White pour l'hétéroscédasticité
Exemple 38 : Diagnostique pour la normalité
Exemple 39 : Construire un modèle
Exemple 40 : Aberrance en X et en Y -- Influence
Télécharger et installer R
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.
Exemple 2 : Affection d'une variable
Exemple 3 : L'éditeur de R
-
Le moyen le plus efficace de travailler en R est d'ouvrir une fenêtre d'édition.
- Sélectionner File -> New Script (en Windows) ou File -> New Document (sur un Mac).
- Nous pouvons écrire un script (une suite de commandes) dans l'éditeur de R.
- Pour soumettre le script, sélectionnez les commandes (avec la souris) que vous souhaitez soumettre.
- Pour soumettre la sélection, appuyez sur CTRL-R (sous Windows) ou appuyez sur CTRL-ENTER (sur un Mac).
- Pour ouvrir une fenêtre d'édition avec Rstudio, sélectionnez File -> New File -> R Script.
- Nous pouvons écrire un script (une séquence de commandes) dans la fenêtre d'édition.
- Pour soumettre le script à la console, sélectionnez les commandes (avec votre souris).
- Pour soumettre la sélection, appuyez sur le bouton Run situé au dessus de la fenêtre d'édition ou appuyez sur CTRL-ENTER (avec Windows) ou appuyez sur CMD-ENTER (avec un Mac).
Exemple 4 : Statistiques Descriptives (univariées)
Soit x un vecteur numérique. Pour notre exemple, considérons l'attribution suivante :
x=c(16,15,14,16,13,12,14,13,10)
- La fonction summary dans R produit un sommaire numérique. Son utilisation est
summary(x)
Voici une sortie de R qui démontre l'affection des valeurs au vecteur numérique et l'utilisation de la fonction summary.
> x=c(16,15,14,16,13,12,14,13,10)
> summary(x)
Min. 1st Qu. Median Mean 3rd Qu. Max.
10.00 13.00 14.00 13.67 15.00 16.00
- Il est possible d'obtenir chacune des statistiques descriptives individuellement.
- Voici une liste de commandes pour certaines statistiques descriptives communes :
sum(x) # pour la somme
mean(x) # pour la moyenne
var(x) # pour la variance
sd(x) # pour l'écart type
min(x) # pour le minimum
max(x) # pour le maximum
median(x) # pour la médiane
quantile(x) # pour certain quantiles
length(x) # pour le nombre de composantes dans le vecteur
- Comme exemple, on calcul la moyenne, l'écart type et certains centiles (minimum, premier quartile, médiane, troisième quartile et
maximum).
> mean(x)
[1] 13.66667
> sd(x)
[1] 1.936492
> quantile(x)
0% 25% 50% 75% 100%
10 13 14 15 16
- Plusieurs méthodes non-paramétriques sont basées sur les rangs. On ordonne les n valeurs du plus petit au plus grand et on utilise
les rangs de 1 à n. Pour des égalilés, on utilise le rang moyen. La fonction rank, nous donne les rangs. Voici un exemple.
> x
[1] 16 15 14 16 13 12 14 13 10
> rank(x)
[1] 8.5 7.0 5.5 8.5 3.5 2.0 5.5 3.5 1.0
Parfois on veut réellement ordonner les valeurs. On peut le faire avec la fonction sort. On obtient les statistiques
d'ordre. Voici un exemple. > x
[1] 16 15 14 16 13 12 14 13 10
> sort(x)
[1] 10 12 13 13 14 14 15 16 16
Remarque: La fonction sort utilise un arrangement en ordre croissant par défaut. On peut utiliser utiliser l'argument
decreasing = TRUE pour obtenir un arrangement en ordre décroissant. Voici un exemple.
> sort(x,decreasing=TRUE)
[1] 16 16 15 14 14 13 13 12 10
Exemple 5 : Statistiques Descriptives (bivariées)
Considérons les deux variables numériques suivantes : y=nombres de rhumes en 6 ans et x=la dose journalière de vitamine C (en mg). Voici
l'attribution des valeurs dans R.
y<-c(16,15,14,16,13,12,14,13,10)
x<-c(0,0,0,25,25,25,50,50,50)
- On peut utiliser la commande
plot(x,y)
pour produire un nuage de points pour décrire l'association entre ces deux variable
numériques. En utilisant la commande suivante
> plot(x,y,xlab="Dose de vitamine C (en mg)",ylab="Nombre de rhumes")
nous avons produit le nuage de point ci bas.
- Pour calculer la covariance et la corrélation entre
x
et y
, on peut utiliser les commandes cov
et cor
.
Voici la sortie de R.
> cov(x,y)
[1] -25
> cor(x,y)
[1] -0.5962848
- Par défaut, R calcul la corrélation de Pearson. Mais on peut ajouter un argument dans la fonction pour obtenir la corrélation de
Spearman. Voici un exemple.
> cor(x,y,method="spearman")
[1] -0.5871366
Remarque: La corrélation de Spearman est la corrélation de Pearson entre les rangs de x et les rangs de y.
Exemple 6 : Loi de probabilités
Loi de probabilités | Nom en R | arguments |
béta | beta | shape1, shape2, ncp |
binomiale | binom | size, prob |
Cauchy | cauchy | location, scale |
khi-deux | chisq | df, ncp |
exponentielle | exp | rate |
F | f | df1, df2, ncp |
gamma | gamma | shape, scale |
géométrique | geom | prob |
hypergéometric | hyper | m, n, k |
log-normale | lnorm | meanlog, sdlog |
logistique | logis | location, scale |
binomiale négative | nbinom | size, prob |
normale | norm | mean, sd |
Poisson | pois | lambda |
rang signé | signrank | n |
t de Student | t | df, ncp |
uniforme | unif | min, max |
Weibull | weibull | shape, scale |
Wilcoxon | wilcox | m, n |
Le tableau ci-haut contient certaines lois de probabilités en R. On utilise le nom de la loi avec un préfixe. Voici la liste des préfixes:
- d pour la densité de probabilité;
- p pour la fonction de répartition;
- q pour la fonction de quantile (l'inverse de la fonction de répartition);
- r pour générer des valeurs au hasard de cette loi.
Voici quelques exemples de son utilisation.
Supposons qu'on veut calculer P(Z>1.65), P(U>12.6), P(T>1.71) et P(F>1.85), où Z suit une loi N(0,1), U suit une loi khi-deux avec 6 degrés de liberté, T suit une loi
t(24) et F suit une loi F(20,15). Voici les commandes.
> 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
Alors, P(Z>1.65)=0.049471, P(U>12.6)=0.0498469, P(T>1.71)=0.05008269, P(F>1.85)=0.1140452.
Le quantile d'ordre p (pour la variable aléatoire X) est une valeur q tel que F(q)=P(X<=q)=p. [Si X est discrète, alors c'est la plus petite quantité q tel que
P(X<=q)>=p.] Par exemple, considérons t(0,975;43), qui est le quantile d'ordre 0,975 de la loi t avec 43 degrés de liberté. Alors, si
T est une variable aléatoire t(43), alors P[T<=t(0,975;43)]=0,975. La commande pour obtenir t(0,975;43) est ci-bas.
> qt(0.975,43)
[1] 2.016692
Alors t(0,975;43)=2.016692.
Supposons qu'on veut générer 10 valeurs d'une loi normale de moyenne 46 et d'écart type 7. Voici la commande.
> 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
Exemple 7 : Importer des données d'un fichier
Nous allons sauvegarder nos données dans un fichier de format CSV (c'est-à-dire un fichier avec des colonnes délimitées par des virgules).
Commentaires:
Exemple 8 : Nuage de points avec une courbe lisse
Exemple: Dans le fichier cerveaux.csv, nous avons n = 25 espèces de mamifères (les rangés). Nous avons deux variables
pour décrire chacune des espèces : masse (en grammes), et gestation (en jours). Nous importons le fichier et nous affichons quelques rangés avec la fonction head
.
> 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
En supposant que x
et y
sont des vecteurs numériques, on peut utiliser la fonction plot(x,y)
pour construire un nuage de point de y
contre x
.
Avec la commande suivante, nous allons construire un nuage de point de la masse du cerveau à la naissance contre la durée de la période de gestation.
with(cerveau, plot(gestation..en.jours.,masse..grammes.,xlab="Gestation (en jours)",ylab="Masse du cerveau (en grammes)"))
Voici le nuage de points.
Afin de mieux visualier le lien entre les deux variables nous pouvons ajouter une courbe lisse au diagramme.
## 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))
Commentaires:
- 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
.
Voici le nuage de points avec la superposition d'une courbe lisse.
Corrélation de Spearman:
-
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 rS=0.86.
> with(cerveau,cor(gestation..en.jours.,masse..grammes.,method="spearman"))
[1] 0.8628308
N.B.: Pour visualiser la linéarité ou non-linéarité d'une association entre deux variables numériques, nous pouvons aussi superimposer une droite estimée sur le nuage de points. On peut utiliser la fonction lm
pour obtenir la droite des moindres carrés et utiliser la fonction abline
pour ajouter cette droite au nuage de points. Pour ajouter la droite au nuage de points et aussi une légende, on utilise les commandes suivantes:
# 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)
Voici le diagramme résultant.
Exemple 9 : Opérations sur une colonne et transformations d'une variable
Soient x
et y
des variables numérique ayant le même nombre de rangés. Soit a
un scalaire (c-à-d un nombre).
Remarque: La plupart des opérations sur les colonnes sont définies par rangée (c'est-à-dire par composante).
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
Voici quelques fonctions pratiques :
sum(x) # la somme des composantes de x
length(x) # le nombre de composantes
Exemple : Considérons l'échantillon aléatoire suivant : 10, 12, 11, 12, 15, 20, 16. Calculons la variance de l'échantillon en utilisant les opérations sur les colonnes. Ensuite vérifions notre calcul avec la fonction var.
> 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
Rappel: La définition de la variance de l'échantillon est
s2=∑ni=1(xi−x¯)2n−1.
Là considérons la variable y=c(11, 12, 11, 13, 14, 18, 17). Calculons la covariance entre x et y avec les opérations sur les colonnes. Ensuite vérifions nos calculs avec la fonction cov.
> sum((x-mean(x))*(y-mean(y)))/(length(x)-1)
[1] 9.404762
> cov(x,y)
[1] 9.404762
Rappel: La définition de la covariance de l'échantillon est x
et y
est
σ^X,Y=∑ni=1(xi−x¯)(yi−y¯)n−1.
Exemple [Carburant] : Considérons les données dans le fichier evaporation.csv. Ce sont des mesures de la vélocité dans l'air (en cm/s) et du coefficient d'évaporation de carburant (en mm2/s)
brûlant dans un moteur d'impulsion. Utilisons la notation suivante : x=la vélocité et y=l'évaporation du carburant.
On importe les données avec R et on affiche quelques rangés.
> 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
Nous calculons les sommes suivantes avec R:
∑i=1nxi=2220; ∑i=1nyi=9,1; ∑i=1ny2i=9,6722; ∑i=1nx2i=580400; ∑i=1nxiyi=2340,4.
> 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
Exemple 10 : Ajuster un modèle linéaire avec la fonction lm
Exemple [Carburant] : Considérons les données dans le fichier evaporation.csv. Ce sont des mesures de la vélocité dans l'air (en cm/s) et du coefficient d'évaporation de carburant (en mm2/s)
brûlant dans un moteur d'impulsion. Utilisons la notation suivante : x=la vélocité et y=l'évaporation du carburant.
On importe les données avec R et on affiche quelques rangés.
> 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
On peut utiliser la fonction lm
pour ajuster le modèle linéaire pour décrire l'évaporation du carbutant un fonction de la vélocité dans l'air.
> mod<-lm(Evaporation ~ Velocite, data=carburant)
> mod
Call:
lm(formula = Evaporation ~ Velocite, data = carburant)
Coefficients:
(Intercept) Velocite
0.059033 0.003807
Alors, le modèle linéaire estimé est Evaporationˆ=0.059033+0.003807velocité.
Commentaires:
- 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, b0=0.059032967 et b1=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=∑ni=1e2i/(n−2)=0.02512. L'écart type de cette variance est dite l'écart type résiduelle se=MSE−−−−√=0.1584 mm2/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
ℓ=−n2ln(2πSSE/n)−n/2=5,758624.
Il y a 3 degrés de liberté, puisque nous avons 3 paramètres à estimer β0,β1 et σ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 R2=89,5%.
> summary(mod)$r.squared
[1] 0.8945677
On peut afficher l'écart type résiduelle se=MSE−−−−√=0,158482.
> summary(mod)$sigma
[1] 0.158482
Exemple 12 : L'analyse de variance avec un objet lm
Exemple: Considérons les données sont dans le fichier suivant: VC.csv. Sauvegarder le fichier dans votre répertoire de travail (current working directory).
Rappel: C'est 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 importe les données avec R et on affiche quelques rangés.
> 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
Utilisons un modèle linéaire simple pour décrire la vitesse de conduction de la moelle épinière en fonction de la taille de l'individu (en cm). La fonction de la moyenne est
E[Y|X=x]=β0+β1x.
Nous voulons tester pour la signification de la taille de l'individu en confrontant
H0: β1=0 contre Ha: β1≠0.
(Puisque nous avons un modèle simple, ceci est aussi dit un test pour la signification de la régression.) Nous pouvons tester ces hypothèses avec la statistique t∗=b1/s{b1}. Mais, une autre approche est d'utiliser une analyse de variance (ANOVA).
Pour obtenir le tableau d'ANOVA, on peut mettre l'objet lm
dans la fonction anova
.
> 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
Commentaires:
- 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 MSR n'estime pas σ2, mais plut\^ot quelque chose de plus grand.
On peut d\'emontrer que
E{MSR}=E{SSR}=σ2(1+β21sxx).
Ceci veut dire que nous avons des preuves que β1≠0.
- Pour avoir une mesure de la signification des preuves contre H0:β1=0, on calcul nos chances d'avoir observer une statistique F aussi grande que 391,3 en supposant que H0 est vraie. Sous H0, F∗∼F(1,n2), 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∗=(b1s{b1})2=SSR/1SSE/(n−2)=F et t21−α/2;n−2=F(1−alpha;1,n−2).
Exemple 13 : Un test F partiel pour comparer deux modèles
Exemple: Considérons les données sont dans le fichier suivant: VC.csv. Sauvegarder le fichier dans votre répertoire de travail (current working directory).
Rappel: C'est 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 importe les données avec R et on affiche quelques rangés.
> 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
Utilisons un modèle linéaire simple pour décrire la vitesse de conduction de la moelle épinière en fonction de la taille de l'individu (en cm). La fonction de la moyenne est
E[Y|X=x]=β0+β1x.
Nous voulons tester pour la signification de la taille de l'individu en confrontant
H0: β1=0 contre Ha: β1≠0.
(Puisque nous avons un modèle simple, ceci est aussi dit un test pour la signification de la régression.) Nous pouvons tester ces hypothèses avec la statistique t∗=b1/s{b1}.
Nous allons vous présenter une autre approche au test dit le test F partiel, qui est une comparaison de modèles. Si H0:β1=0 est vraie, alors nous imposons une contrainte \`a l'espace des param\`etres et ceci réduit le modèle.
La fonction de la moyenne du modèle réduit est
E[Y|X=x]=β0.
Alors, voici une autre façon d'écrire les hypothèses\:
H0: E[Y|X=x]=β0 contre Ha: E[Y|X=x]=β0+β1x.
Pour évaluer les preuves contre H0 en faveur de Ha, nous allons comparer l'ajustement des deux modèles tel que décrit par la somme des carrés résiduelles. Avec R, on ajuste les deux modèles et pour les comparer on met les deux objets dans la fonction anova
.
> 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
Nous comparons la somme des carrés du modèle complet SSE=112,44 à la somme des carrés du modèle réduit SSE(R)=400,00, en calculant la somme des carrés supplémentaires
ExtraSS=SSE(R)−SSE=287.56.
Plus que cette somme de carrés supplémentaire est grande, plus que nous considérerons le modèle réduit comme étant moins bien ajusté en comparaison au modèle complet. Si ces preuves sont fortes, alors nous avons des preuves contre H0 en faveur de Ha que la pente est non-nulle.
Exemple 14 : Une variable explicative binaire
Exemple [Azote] : Une étude fut menée sur le développement de l'ectomycorhizue, une relation symbiotique entre
les racines des arbres et un champignon, dans lequel les minéraux sont transférés du champignon aux arbres et du
sucre des arbres au champignon, 20 semis de chênes rouges du nord exposés au champignon pisolithus tinctorus ont
été cultivés dans une serre. Toutes les semis ont été plantées dans le même type de terre et ont reçu la même quantité
de soleil et d'eau. La moitié des semis (choisis au hasard) ont reçu un traitement de 368 ppm d'azote sous la forme
NaNO3. Les autres semis n'ont pas reçu d'azote. La masse de la tige, en grammes, est mesurée après 140 jours. Les
données sont dans le fichier Azote.csv.
Nous importons les données et nous affichons quelques rangés. On remarque que la variable Azote
est une variable catégorique.
> 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
Affichons les catégories de la variable Azote
. R considère la première catégorie comme la catégorie de référence. Alors, ne pas avoir reçu de l'azote est considéré comme le groupe de référence. On peut changer l'ordre des catégories afin que le groupe de traitement d'azote soit le groupe de référence avec la fonction factor
.
> levels(azote$Azote)
[1] "non" "oui"
> azote$Azote<-factor(azote$Azote, levels=c("oui","non"))
> levels(azote$Azote)
[1] "oui" "non"
Commentaire: Si nous avons un groupe de contrôle, c'est ce groupe qui est souvent utilisé comme le groupe de référence. Alors, nous allons changer l'ordre afin que le groupe sans azote soit le groupe le premier groupe.
> azote$Azote<-factor(azote$Azote, levels=c("non","oui"))
Voici des statistiques descriptives pour la masse de la tige (y) pour chacun des groupes.
> 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
En alternatif, on peut utiliser la fonction tapply
.
> 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
Remarque: Avec la commande tapply(y,x,FUN)
, on évalue la fonction FUN
sur la variable y
selon les niveaux de la variable x
.
En supposant que le groupe de contrôle soit le groupe 1 et que le groupe de traitement soit le groupe 2, on a
y¯1=0,392; y¯2=0,628; s21=0,013551; s22=0,018196; n1=n2=10.
Si on suppose que les variances des deux populations sont égales, alors on peut estimer cette variance avec la variance pondérée
s2p=(n1−1)s21+(n2−1)s22n1+n2−2=0,014286.
Avec R :
> ## variance poindérée
> sum(variance*(n-1)/sum(n))
[1] 0.014286
Pour comparer les deux groupes visuellement, on peut utiliser des diagrammes à boîte comparatifs avec une superposition de points.
boxplot(Masse~Azote,data=azote,ylab="Masse",
names=c("Sans azote","Avec azote"))
stripchart(Masse~Azote,data=azote,method="jitter",
vertical = TRUE,add=TRUE)
Voici le diagramme correspondant.
On peut utiliser un test T de Student pour comparer les deux moyennes. La valeur observée de la statistique du test est
t∗=y¯1−y¯2sp1/n1+1/n2−−−−−−−−−−√=−4,1885.
La valeur p est 2P(T(18)>|−4,1885|)=0,0005521.
Avec R :
> 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
On peut utiliser une approche de régression pour faire ce test. (Ceci va nous aider à généraliser le test à la comparaison de plus que 2 groupes plus tard.)
Pour identifier les groupes, nous allons utiliser une variable muette (dummy variable en anglais)~:
xi={1,0,group de traitementsinon.
Consid\'erons le modèle de régression linéaire simple: Y1,Y2,…,Yn sont normales, ind\'ependantes tel que
E[Yi]=β0+β1xi={β0=μ1,β0+β1=μ2,xi=0xi=1
et $V[Y_i]=\sigma^2. Alors, nous avons deux populations normales indépendantes avec des variances égales.
Quand notre variable explicative est une variable catégorique (ceci est dit factor dans la terminologie de R), alors R va coder des variables muettes correspondantes. Pour afficher ces variables muettes, on utilise la fonction contrasts
.
> contrasts(azote$Azote)
oui
non 0
oui 1
On interprète la variable de la façon suivante. Ceci est utilisé pour identifier la catégorie oui
. Elle est égale à 0, si l'unité est non. Elle est égale à 1, si l'unité est oui.
Ajustons le modèle linéaire pour décrire la masse selon le groupe de traitement.
> mod<-lm(Masse~Azote,data=azote)
> mod$coefficients
(Intercept) Azoteoui
0.392 0.236
> summary(mod)$sigma^2
[1] 0.01587333
Alors, notre estimation de la masse moyenne est
μ^Y|X=x=b0+b1x={b0=0.392,b0+b1=0.392+0.236=0.628,x=0x=1
L'estimation de la variance σ2 est MSE=0.01587333. Remarquons que ceci nous donne y¯1, y¯2, et s2p, respectivement.
Remarquons que μ1=μ2 si et seulement si β1=0. Alors tester pour la signification de la variable explicative qui identifie le traitement est interpréter comme un test pour l'égalité des moyennes. La valeur observée de la statistique du test est
t∗=b1/s{b1}=4,1885 et la valeur p du test est 2P(T(18)>|4,1885|)=0,0005521.
> 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
Exemple 15 : Inférence concernant une corrélation
Exemple [VC] : Considérons les données sont dans le fichier suivant: VC.csv.
Sauvegarder le fichier dans votre répertoire de travail (current working directory).
Rappel: C'est 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 importe les données avec R et on affiche quelques rangés.
> 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
Supposons que nous voulons testez H0:ρ=0 (o\`u ρ est la corrélation (de Pearson) entre la vitesse de conduction et la taille de l'individu) contre Ha:ρ≠0. On peut utiliser la fonction cor.test
.
> 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
Alors, la valeur observée de la statistique du test est
t∗=rn−21−r2−−−−−−√=19,781.
La valeur p est 2P(T(153)>|19,781|)<0,00001. Les preuves que la corrélation est non-nulle sont significatives.
Commentaires:
- Pour obtenir l'intervalle de confiance pour la corrélation, la fonction
cor.test
utilise la transformation z de Fisher. Cette transformation est
12ln(1+r1−r)=arctanh(r).
Fisher a démontré que la distribution d'échantillonage de cette transformation suit approximativement
N(12ln(1+ρ1−ρ),1n−3).
Alors, un IC \`a 95\% pour arctanh(ρ) est
12ln(1+x1−x)±z0,975n−3−−−−−√=[a,b],
où z0,975=1,96. Ensuite en appliquant la fonction inverse qui est
tanh(x)=e2x−1e2x+1,
on obtient un IC pour ρ. Alors, un IC \`a 95\% pour ρ est
[e2a−1e2a+1,e2b−1e2b+1]=[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
Exemple 21 : Introduction au modèle de régression linéaire multiple (aussi dit modèle linéaire)
Exemple [Examen SAT en 1999]: Considérons les données de une étude des résultats aux examens SAT aux Etats-Unis en 1999.
Les données sont dans le fichier : sat.csv. On importe les données et on affiche quelques rangés du jeu de données.
> 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
Commentaires: Il y a 50 unités statistiques (les rangés). Les unités sont les 50 états aux Etats-Unis. Nous avons 8 variables pour décrire ces unités. Les trois premières variables sont 3 variables définie au niveau de l'état qui pourrait
être lié aux résultats aux examens du SAT.
Nous voulons étudier le lien entre le revenu moyen des enseignents et les résultats aux SAT (test standardisé
utilisé pour les admissions au collège aux Etats-Unis) en 1999. On observe que l'association entre le salaire moyen et
les résultats aux SAT est significative (p = 0,0031). Cette association est approximativement linéaire, négative avec
une corrélation de r=−R2−−−√=−0,44.
Avec R :
> 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
Question : Est-ce que nous avons identifié un lien causal? Est-ce que le salaire élevé des enseignants de certains
états est la cause des résultats médiocres aux tests? Est-ce que certains états devraient réduire le salaire des enseignants
pour améliorer la performance des étudiants aux tests?
Réponse : Non, quand nous observons une association significative entre deux variables, ceci ne veut pas
nécessairement impliquer qu'une est la cause de l'autre. Il peut exister une autre variable qui est liée aux deux
variables et qui explique l'assocation qu'on ait observée. Si une telle variable existe, elle est dite une variable
confusionnelle (confounding variable en anglais).
Variable confusionnelle possible : Essayons de penser à une explication du phénomène ci-haut qui va contre
le bon sense. Il est possible que le taux de participation au test soit lié aux résultats et aux salaires.
Stratification des unités avec une variable muette : La médiane du taux de participation est 28%, alors
nous allons utiliser ce niveau pour définir une variable muette :
MuetteHaut={1;0;taux de participation ≥28%taux de participation <28%
Nous avons diviser les états en 2 groupes : haut taux de participation et faible taux de participation.
Avec R :
> 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)
Commentaire:
- 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))
Un modèle linéaire pour estimer les deux droites: Considérons un modèle linéaire avec 3 variables
indépendantes :
x1=MuetteHaut; x2=Salaire moyen; x3=MuetteHaut * Salaire moyen.
La fonction de la moyenne pour ce mod\`ele est
E{Y}==β0+β1x1+β2x2+β3x3{β0+β1+(β2+β3)(salaire moyen),β0+β2(salaire moyen),MuetteHaut =1MuetteHaut =0
L'estimation par les moindres carrés nous donne:
Strate |
Ordonnée à l'origine |
Pente |
MuetteHaut = 1 |
b0+b1=861,05 |
b2+b3=1,0657 |
MuetteHaut = 0 |
b0=1035,71 |
b2=0,1733 |
Avec R:
> ## 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
Revoici le nuage de points avec la superposition des deux droites.
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)
Taux de participation comme variable quantitative : Dans la pratique, si la variable confusionnelle est une
variable quantitative, alors on devrait l'inclure dans le modèle comme une des variables indépendantes. C'est-à-dire,
on ne devrait pas transformer une variable quantitative en variable catégorique. La plupart du temps ceci cause une dimimution de précision et de puissance.
Nous utilisons un modèle linéaire ou le résultat total aux tests SAT est la variable dépendante et les variables indépendantes
sont le salaire moyen et le taux de participation. La valeur explicative du modèle est décrite avec R2=80,56% et
R2a=79,73%. L'écart type résiduel est se=33,69. L'association entre le taux de participation et les résultats aux tests SAT
est significative et négative. Tandis que l'association entre le salaire moyen des enseignants et le résultat au SAT est
significative et positive.
> 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
Le modèle linéaire estimée est
Total SATˆ=987,9005+2,1804(Salaire Moyen)−2,7787(Taux de participation (en %)).
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×p. Supposons que nous voulons tester H0 qui nous donne un modèle réduit qui est emboité dans notre modèle. Supposons que X0 de taille n×q est la matrice du plan du
modèle réduit, où les colonnes de X0 appartiennent à l'espace colonne de X. On peut tester H0 avec la statistique F suivante:
F∗=ExtraSS/(p−q)SSE/(n−p), où ExtraSS=SSE(R)−SSE.
Plus que ExtraSS est grande, plus que nous avons des preuves contre H0. Sous H0, 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ée 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"
Remarques:
Nous allons utliser un modèle linéaire pour décrire le temps de rétablissement selon la forme. La fonction de la moyenne du modèle linéaire est
\[
E\{Y\}=\beta_0+\beta_1\,\,I\{\mbox{Groupe}=2\} +\beta_2\,I\{\mbox{Groupe}=3\}.
\]
Ici on utilise la notation \(I\) pour une fonction indicatrice, où \(I(A)=1\), si \(A\) est réalisé, sinon \(I(A)=0\). Alors, on a
\[
\begin{eqnarray*}
E\{Y\} &=&\beta_0+\beta_1\,\,I\{\mbox{Groupe}=2\} +\beta_2\,I\{\mbox{Groupe}=3\} \\
&=& \left\{
\begin{array}{ll}
\beta_0=\mu_1, & \mbox{Groupe}=1\\
\beta_0+\beta_1=\mu_2 & \mbox{Groupe}=2\\
\beta_0+\beta_2=\mu_3 & \mbox{Groupe}=3\\
\end{array}
\right.
\end{eqnarray*}
\]
Remarques:
- 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.
Nous ajustons le modèle d'ANOVA pour décrire le temps de rétablissement selon la forme du patient. Le modèle est significatif \((F(2,21)=16.96; p<0.0001)\).
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
Commentaires:
- 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}=β0+β1x1+…+βp−1xp−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
H0 : βj=0, pour j=1,2,…,p−1 contre Ha : βj≠0, pour au moins un j .
La statistique du test est
F∗=MSRMSE=SSR/(p−1)SSE/(n−p),
qui suit une loi F(p−1,n−p) si H0 est vraie. Des grandes valeurs de F∗ sont considéré comme des preuves contre H0.
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 x1, et
le revenu disponible par habitant x2 (en milliers de dollars).
Le mod\`ele de r\'egression estim\'ee est
y^=−68,8571+0,00145x1+9,3655x2.
> 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 x1 et x2 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 H0:β1=β2=0 contre Ha: au moins un des coefficients β1,β2 n'est pas zéro.
La statistique du test est
F∗=SSR/(p−1)SSE/(n−p)=99,1
et la valeur-p est P(F(2,18)>99,1)=1,921×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 H0:β1=β2=0. Ceci est équivalent à confronter
H0:E{Y}=β0 contre Ha:E{Y}=β0+β1x1+β2x2.
Nous ajustons les deux modèles et on les compare avec la fonction anova
. La statistique du test est
F∗=Extra/(p−q)SSE/(n−p)=99,103,
et la valeur-p est P(F(2,18)>99,103)=1,92×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,λ),
o\`u
λ=E{Y}′(H−H0)E{Y}σ2=nf2,
et H est la matrice chapeau du mod\`ele complet et H0 est la matrice chapeau du mod\`ele r\'eduit, et f2 est le f2 de Cohen.
Voici un graphe de la densité de probabilité de la loi F(p−q=5,n−p=13,λ), pour λ=0, λ=1,\lambda=2,et\lambda=3.NousavonsaussisuperimposéunedroiteverticaleàlavaleurcritiqueF(0,\!95;5,13).Rappel:OnrejetteH_0siF^*>F(0,\!95;5,13).Plusque\lambda=n\,f^2estgrand,plusquechancederejetterH_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 H0. Quand H0 est vraie, la puissance est P(rejet H0|H0 est vraie)=α. Mais, si Ha 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é λ=nf2, où f2 est le f2 de Cohen. La puissance du test linéaire générale est
puissance=P(F(p−q,n−p,λ=nf2)>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) f2=0,1; (b) f2=2,5.
(a) La puissance est
puissance=P(F(p−q=5,n−p=13,λ=nf2=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
puissance=P(F(p−q=5,n−p=13,λ=nf2=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 f2=1,5 avec une puissance de 80\% \`a un niveau de signification de α=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 f2=2,5, est
puissance=P(F(p−q=5,n−p=15−7=8,λ=nf2=(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×p. Supposons que nous voulons tester H0 qui nous donne un modèle réduit qui est emboité dans notre modèle. Supposons que X0 de taille n×q est la matrice du plan du
modèle réduit, où les colonnes de X0 appartiennent à l'espace colonne de X. On peut tester H0 avec la statistique F suivante:
F∗=ExtraSS/(p−q)SSE/(n−p), où ExtraSS=SSE(R)−SSE.
Plus que ExtraSS est grande, plus que nous avons des preuves contre H0. Sous H0, 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}=β0+β1I{Groupe=2}+β2I{Groupe=3}+β3(Age)+β4(Age)I{Groupe=2}+β5(Age)I{Groupe=3}.
Nous voulons tester l'hypothèse nulle qu'il n'y a pas d'interactions entre la forme et l'âge.
H0 : E{Y}=β0+β1I{Groupe=2}+β2I{Groupe=3}+β3(Age)
contre
Ha : E{Y}=β0+β1I{Groupe=2}+β2I{Groupe=3}+β3(Age)+β4(Age)I{Groupe=2}+β5(Age)I{Groupe=3}.
On ajuste les deux modèles et on les compare avec la fonction anova
. N.B. En imposant les contraintes de H0, 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∗=ExtraSS/(p−q)SSE/(n−p)=0,22184/25,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}=β0+β1I{Groupe=2}+β2I{Groupe=3}+β3(Age).
- Test pour la signification de l'âge: On veut tester H0:β3=0 contre Ha:β3≠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∗=b3/s{b3}=36,4608. La valeur-p est 2P(T(20)>36,4608)=9,1×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
H0 : E{Y}=β0+β1I{Groupe=2}+β2I{Groupe=3} contre Ha : E{Y}=β0+β1I{Groupe=2}+β2I{Groupe=3}+β3(Age).
La statistique du test est
F∗=ExtraSS/(p−q)SSE/(n−p)=1329,4.
La valeur-p du test est P(F(1,20)>1329,4)=2,2×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}=β0+β1I{Groupe=2}+β2I{Groupe=3}+β3(Age).
On veut tester
H0:β1=0,β2=0 contre Ha:H0 est fausse.
Ici on ne peut pas utiliser un test t. Mais, on peut utiliser un test F. On veut tester
H0 : E{Y}=β0+β3(Age) contre Ha : E{Y}=β0+β1I{Groupe=2}+β2I{Groupe=3}+β3(Age).
La statistique du test est
F∗=ExtraSS/(p−q)SSE/(n−p)=399,11.
La valeur-p du test est P(F(2,20)>399,11)=2,2×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
H0 : E{Y}=β0 contre Ha : E{Y}=β0+β1I{Groupe=2}+β2I{Groupe=3}.
et les hypothèses pour Age sont
H0 : E{Y}=β0+β1I{Groupe=2}+β2I{Groupe=3}. contre Ha : E{Y}=β0+β1I{Groupe=2}+β2I{Groupe=3}+β3(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
H0 : E{Y}=β0. contre Ha : E{Y}=β0+β1(Age).
et les hypothèses pour Groupe sont
H0 : E{Y}=β0+β1(Age). contre Ha : E{Y}=β0+β1(Age)+β2I{Groupe=2}+β3I{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}=β0+β1x1+…+βp−1xp−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
H0 : βj=0, pour j=1,2,…,p−1 contre Ha : βj≠0, pour au moins un j .
La statistique du test est
F∗=MSRMSE=SSR/(p−1)SSE/(n−p),
qui suit une loi F(p−1,n−p) si H0 est vraie. Des grandes valeurs de F∗ sont considéré comme des preuves contre H0.
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 x1, et
le revenu disponible par habitant x2 (en milliers de dollars).
Le mod\`ele de r\'egression estim\'ee est
y^=−68,8571+0,00145x1+9,3655x2.
> 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 x1 et x2 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 H0:β1=β2=0 contre Ha: au moins un des coefficients β1,β2 n'est pas zéro.
La statistique du test est
F∗=SSR/(p−1)SSE/(n−p)=99,1
et la valeur-p est P(F(2,18)>99,1)=1,921×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 H0:β1=β2=0. Ceci est équivalent à confronter
H0:E{Y}=β0 contre Ha:E{Y}=β0+β1x1+β2x2.
Nous ajustons les deux modèles et on les compare avec la fonction anova
. La statistique du test est
F∗=Extra/(p−q)SSE/(n−p)=99,103,
et la valeur-p est P(F(2,18)>99,103)=1,92×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,λ),
o\`u
λ=E{Y}′(H−H0)E{Y}σ2=nf2,
et H est la matrice chapeau du mod\`ele complet et H0 est la matrice chapeau du mod\`ele r\'eduit, et f2 est le f2 de Cohen.
Voici un graphe de la densité de probabilité de la loi F(p−q=5,n−p=13,λ), pour λ=0, λ=1,\lambda=2,et\lambda=3.NousavonsaussisuperimposéunedroiteverticaleàlavaleurcritiqueF(0,\!95;5,13).Rappel:OnrejetteH_0siF^*>F(0,\!95;5,13).Plusque\lambda=n\,f^2estgrand,plusquechancederejetterH_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 H0. Quand H0 est vraie, la puissance est P(rejet H0|H0 est vraie)=α. Mais, si Ha 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