Ordre d'introduction des effets principaux et interactions

[François Bonnot]

Bonjour,
Il y a eu récemment une question sur la syntaxe aov :
viewtopic.php?t=587
qui me conduit à la question suivante.

Supposons un essai avec 2 facteurs de contrôle control1 et control2, 1 facteur étudié etud, et une variable x. Le jeu de données fictif suivant en est une illustration dans un cas volontairement déséquilibré :

Code : Tout sélectionner

df <- data.frame(control1=as.factor(c(1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2)),
                 control2=as.factor(c(1, 1, 1, 1, 2, 2, 1, 1, 1, 2, 2, 2)),
                 etud=as.factor(c(1, 2, 3, 3, 1, 2, 2, 2, 3, 1, 3, 3)),
                 x=c(10.22, 11.16, 11.04, 10.72,  3.21,  4.78,  3.80,
                   5.39,  5.31,  4.27,  4.68,  4.76))
df$control1.2 <- as.factor(paste(as.numeric(df$control1),as.numeric(df$control2),sep='.'))


Le facteur control1.2 n'est autre que la juxtaposition de control1 et control2, et correspond donc à un facteur de contrôle équivalent aux 2 premiers.
L'analyse sans tenir compte des facteurs de contrôle ne montre pas d'effet de etud :

Code : Tout sélectionner

> anova(lm(x~etud,data=df))
Analysis of Variance Table

Response: x
          Df  Sum Sq Mean Sq F value Pr(>F)
etud       2   4.336   2.168  0.1867 0.8329
Residuals  9 104.522  11.614               


Si on fait rentrer les facteurs de contrôle par le biais de control1.2, on n'obtient pas non plus d'effet de etud, malgré une diminution drastique de la résiduelle :

Code : Tout sélectionner

> anova(lm(x~control1.2+etud,data=df))
Analysis of Variance Table

Response: x
           Df  Sum Sq Mean Sq  F value    Pr(>F)   
control1.2  3 105.353  35.118 101.1405 1.583e-05 ***
etud        2   1.421   0.711   2.0467    0.2101   
Residuals   6   2.083   0.347                       


Maintenant on veut intégrer explicitement control1 et control2 dans le modèle, et c'est là que les problèmes commencent :

Code : Tout sélectionner

> anova(lm(x~control1*control2+etud,data=df))
Analysis of Variance Table

Response: x
                  Df Sum Sq Mean Sq F value    Pr(>F)   
control1           1 43.777  43.777 126.080 2.981e-05 ***
control2           1 31.507  31.507  90.743 7.635e-05 ***
etud               2  8.454   4.227  12.174 0.0077288 **
control1:control2  1 23.036  23.036  66.346 0.0001841 ***
Residuals          6  2.083   0.347                     


En effet l'interaction control1:control2 est rejetée en fin d'analyse (et est effectivement introduite à la fin dans les calculs), de sorte que le facteur etud est significatif à tort car il absorbe une partie de l'interaction des facteurs de contrôle.
L'analyse correcte serait la suivante :

Code : Tout sélectionner

> anova(lm(x~control1+control2+control1.2+etud,data=df))
Analysis of Variance Table

Response: x
           Df Sum Sq Mean Sq  F value    Pr(>F)   
control1    1 43.777  43.777 126.0799 2.981e-05 ***
control2    1 31.507  31.507  90.7425 7.635e-05 ***
control1.2  1 30.069  30.069  86.5991 8.716e-05 ***
etud        2  1.421   0.711   2.0467    0.2101   
Residuals   6  2.083   0.347                       


D'où la question : est-il possible d'imposer l'odre d'introduction des effets et interactions dans un modèle ? par défaut il semble que R rejette systématiquement les intéractions à la fin.

[Nicolas Péru]

En fait je me demande si la question posée ne serait pas un peu plus pertinente s'il y avait un effet de la variable etud? A cause de cette absence d'effet, je me demande si dans le dernier modèle (étude placée en dernier on ne vient pas juste vérifier la possibilité d'une interaction control1xetud ou control2xetud ce qui serait plus explicite avec le modèle complet?
A priori, on ne peut pas vraiment imposer autrement l'ordre des effets qu'en bricolant une interaction. Ceci vient peut être de ce que l'on fait en général de l'interaction...on espère franchement qu'elle ne soit pas significative :)

La question de l'ordre des variable ne m'a jamais paru très claire mais bon...je me suis contenté jusque là de l'explication du plan deséquilibré sans pouvoir (chercher à:D) la cerner vraiment.

[Tillard]

Bonjour
une question similaire a déjà été débattue sur le forum de R

voir http://finzi.psych.upenn.edu/R/Rhelp02a ... 14240.html

Code : Tout sélectionner

> anova(aov(x ~ control1 * control2 + etud, data=df))
Analysis of Variance Table

Response: x
                  Df Sum Sq Mean Sq F value    Pr(>F)   
control1           1 43.777  43.777 126.080 2.981e-05 ***
control2           1 31.507  31.507  90.743 7.635e-05 ***
etud               2  8.454   4.227  12.174 0.0077288 **
control1:control2  1 23.036  23.036  66.346 0.0001841 ***
Residuals          6  2.083   0.347                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> anova(aov(terms(x ~ control1 * control2 + etud, keep.order=T),df))
Analysis of Variance Table

Response: x
                  Df Sum Sq Mean Sq  F value    Pr(>F)   
control1           1 43.777  43.777 126.0799 2.981e-05 ***
control2           1 31.507  31.507  90.7425 7.635e-05 ***
control1:control2  1 30.069  30.069  86.5991 8.716e-05 ***
etud               2  1.421   0.711   2.0467    0.2101   
Residuals          6  2.083   0.347                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cordialement

[François Bonnot]

En fait je me demande si la question posée ne serait pas un peu plus pertinente s'il y avait un effet de la variable etud?

Non, au contraire : le fait qu'etud soit (faussement) significatif dans un cas et pas dans l'autre accentue même le problème. Mais l'analyse sortie par défaut est incorrecte, que l'effet de etud soit ou non significatif.

A cause de cette absence d'effet, je me demande si dans le dernier modèle (étude placée en dernier on ne vient pas juste vérifier la possibilité d'une interaction control1xetud ou control2xetud ce qui serait plus explicite avec le modèle complet?

Non : on place volontairement etud en dernier pour tester son effet conditionnellement à tous les degrés de liberté des facteurs de contrôle (donc y compris leur interaction). Si on ne fait pas cela, le déséquilibre de l'essai fait que la somme des carrés calculée pour etud (donc placée avant l'interaction) provient à la fois de etud et d'une partie de l'interaction control1:control2 (pour s'en persuader on peut vérifier que la somme des carrés etud+control1:control2 est constante bien que décomposée différemment dans les 2 cas : 8.454+23.036=30.069+1.421=31.490). On ne teste alors pas etud, mais un mélange de etud et de l'interaction, c'est pour cela que l'analyse est incorrecte.

[François Bonnot]

Merci Emmanuel pour ce lien.
Si keep.order n'avait pas existé, il aurait fallu l'inventer.

[Logez Maxime]

Bonjour,


Je crois qu'ici on est en droit de se poser quelques questions tout de même :
quel importance attribué à une analyse de variance sur 12 valeurs avec 3 facteurs qui prennent respectivement 2 ou 3 modalités, et plus que 6 df pour estimer la variance résiduelle ?
Quelle signification ce SCEres peut-il avoir ?
De plus si etud est pas significatif au début pourquoi vouloir le conserver par la suite ? En général on a plutôt tendance à utilisé le principe de parcimonie s'il n'est pas significatif on l'enlève et on travail sur les deux autres facteurs et après on regarde de possibles interactions.

Maxime

[François Bonnot]

Je crois qu'ici on est en droit de se poser quelques questions tout de même :
quel importance attribué à une analyse de variance sur 12 valeurs avec 3 facteurs qui prennent respectivement 2 ou 3 modalités, et plus que 6 df pour estimer la variance résiduelle ?
Quelle signification ce SCEres peut-il avoir ?

Le jeu de données n'a pas de valeur statistique, il est fictif et je l'ai créé en y injectant les effets que je voulais (ou non) voir apparaître pour illustrer la question sur l'ordre des effets. Il est de petite taille mais suffisant sur le plan numérique puisqu'il permet de vérifier que l'analyse proposée par Emmanuel est exactement la même que celle que je voulais obtenir. Evidemment dans la vraie vie on a des jeux de données plus étoffés, mais le principe reste le même.

De plus si etud est pas significatif au début pourquoi vouloir le conserver par la suite ?

En principe il faut définir du premier coup un modèle d'analyse cohérent avec le plan d'expérience réalisé et la question scientifique que l'on se pose. Ici le bon modèle est
x ~ control1 + control2 + control1.2 + etud, que l'on écrira donc plutôt
x ~ control1 * control2 + etud, keep.order=T pour éviter d'avoir à former le facteur control1.2

Les autres modèles ont été présentés pour illustrer le problème technique, mais en principe il ne faut pas les utiliser. En particulier le modèle x~etud n'a aucun intérêt (je l'ai fait figurer pour mettre en évidence la réduction de l'erreur après prise en compte des facteurs de contrôle), il est même dangereux car il peut conduire à des conclusions fausses dans les 2 sens (erreur de 1ère ou de 2ème espèce)

En général on a plutôt tendance à utilisé le principe de parcimonie s'il n'est pas significatif on l'enlève et on travail sur les deux autres facteurs et après on regarde de possibles interactions.

Dans ce cas particulier cela ne présente pas d'intérêt puisque etud est le seul facteur étudié (donc l'objet même de la question scientifique que l'on se pose), on ne va tout de même pas analyser les facteurs de contrôle à vide! (n'oublions pas qu'ils ne sont là que pour diminuer la résiduelle)

Dans un cas plus général on n'a pas intérêt à procéder comme ça, car si les F des effets non significatifs sont supérieurs à 1 (ce qui est souvent le cas) on risque, en les retirant, de perdre de la puissance sur les autres effets, sauf si l'erreur résiduelle est estimée avec un petit nombre de ddl.

[Nicolas Péru]

merci pour ces précisions sur mon post François :)

Mais ceci dit la question de la puissance du test est tout de même importante dans ces cas là non ?
n= 12 pour tester l'effet de 3 facteurs à respectivement 2,2 et 3 niveaux ? Sans parler des interactions...
Cela ne peut il pas affecter les calculs de p-value ? Peut on réellement considérer chacune de ces analyses comme valable ?

Mais effectivement, peut être que cela ne change pas grand chose fondamentalement sur le principe de l'ordre d'introduction des variables.

petite poursuite du sujet :) :

que penser de ça alors ?

Code : Tout sélectionner

>summary(glm(x~control1*control2+etud,data=df))

Call:
glm(formula = x ~ control1 * control2 + etud, data = df)

Deviance Residuals:
    Min       1Q   Median       3Q      Max 
-1.0411  -0.1985   0.0637   0.2878   0.5489 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)   
(Intercept)           10.140      0.435   23.32  4.1e-07 ***
control12             -6.174      0.485  -12.72  1.4e-05 ***
control22             -6.583      0.546  -12.07  2.0e-05 ***
etud2                  0.875      0.510    1.71  0.13723   
etud3                  0.852      0.463    1.84  0.11511   
control12:control22    6.619      0.813    8.15  0.00018 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 0.3472)

    Null deviance: 108.8580  on 11  degrees of freedom
Residual deviance:   2.0833  on  6  degrees of freedom
AIC: 27.04

Number of Fisher Scoring iterations: 2

#################################################

> summary(glm(x~control1+control2+control1.2+etud,data=df))

Call:
glm(formula = x ~ control1 + control2 + control1.2 + etud, data = df)

Deviance Residuals:
    Min       1Q   Median       3Q      Max 
-1.0411  -0.1985   0.0637   0.2878   0.5489 

Coefficients: (2 not defined because of singularities)
              Estimate Std. Error t value Pr(>|t|)   
(Intercept)    10.1404     0.4349   23.32  4.1e-07 ***
control12      -6.1743     0.4854  -12.72  1.4e-05 ***
control22       0.0361     0.5596    0.06  0.95070   
control1.21.2  -6.6189     0.8126   -8.15  0.00018 ***
control1.22.1       NA         NA      NA       NA   
control1.22.2       NA         NA      NA       NA   
etud2           0.8750     0.5103    1.71  0.13723   
etud3           0.8518     0.4625    1.84  0.11511   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 0.3472)

    Null deviance: 108.8580  on 11  degrees of freedom
Residual deviance:   2.0833  on  6  degrees of freedom
AIC: 27.04

Number of Fisher Scoring iterations: 2


Vu ces résultats, et le fait que l'option keep.order n'est pas dispo dans anova(glm()), cela veut il dire que la manière dont procède R dans l'analyse de variance d'un glm, conserve l'ordre comme si keep.order était par défaut laissé à TRUE ?

Il me semble pourtant bien qu'un glm avec family=gaussian (valeur par défaut) revient à un lm donc pourquoi cette différence ???...pas très clair tout ça. :\

[Tillard]

Bonjour

en reprenant les données fictives de Francois, je n'obtiens pas la même chose avec

Code : Tout sélectionner

> model2 <- glm(terms(x ~ control1 * control2 + etud, keep.order=T), data=df)
> summary(model2)

Call:
glm(formula = terms(x ~ control1 * control2 + etud, keep.order = T),
    data = df)

Deviance Residuals:
     Min        1Q    Median        3Q       Max 
-1.04107  -0.19848   0.06375   0.28777   0.54893 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)   
(Intercept)          10.1404     0.4349  23.318 4.08e-07 ***
control12            -6.1743     0.4854 -12.720 1.45e-05 ***
control22            -6.5829     0.5455 -12.067 1.97e-05 ***
control12:control22   6.6189     0.8126   8.145 0.000184 ***
etud2                 0.8750     0.5103   1.715 0.137234   
etud3                 0.8518     0.4625   1.842 0.115110   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 0.3472179)

    Null deviance: 108.8580  on 11  degrees of freedom
Residual deviance:   2.0833  on  6  degrees of freedom
AIC: 27.043

Number of Fisher Scoring iterations: 2


je ne vois pas effectivement d'ou peut venir la difference entre ces résultats et ceux que vous obtenez en utilisant ?

Code : Tout sélectionner

> summary(glm(x~control1+control2+control1.2+etud,data=df))

???????????

cordialement

[Nicolas Péru]

je ne vois pas effectivement d'ou peut venir la difference entre ces résultats et ceux que vous obtenez en utilisant ?
Code:

> summary(glm(x~control1+control2+control1.2+etud,data=df))



???????????

justement je crois qu'il n'y en a pas :)
C'est justement ce qu m'étonne un peu. Ou plutôt ce que je n'arrive pas bien à comprendre. Votre model2 correspond à mon premier modèle sauf que je n'emploie pas la fonction terms() et que du coup l'introduction de keep.order=T me renvois un message d'erreur du type terme inutilisé.
Par contre si on regarde vos valeurs de paramètres et les miennes ou l'AIC ou quoique ce soit c'est rigoureusement identique à l'arrondi près.

Dans mon second modèle, R gère l'interaction introduite comme un nouveau factueur et va donc tester les interactions entre les niveaux de ce facteur...ce n'est pas vraiment ce que l'on veut obtenir je pense...

[François Bonnot]

Mais ceci dit la question de la puissance du test est tout de même importante dans ces cas là non ?
n= 12 pour tester l'effet de 3 facteurs à respectivement 2,2 et 3 niveaux ? Sans parler des interactions...
Cela ne peut il pas affecter les calculs de p-value ? Peut on réellement considérer chacune de ces analyses comme valable ?

Bonjour,
Effectivement, il faudrait prendre les précautions d'usage sur la puissance pour interpréter ces résultat s'il sagissait d'un vrai essai...

...le fait que l'option keep.order n'est pas dispo dans anova(glm())...

En fait c'est la fonction term(), plus précisément terms.formula, qui gère keep.order (voir le dernier message d'Emmanuel), donc ça doit marcher avec n'importe quelle fonction stat qui utilise les modèles.

Vu ces résultats ... cela veut il dire que la manière dont procède R dans l'analyse de variance d'un glm, conserve l'ordre comme si keep.order était par défaut laissé à TRUE ?

en reprenant les données fictives de Francois, je n'obtiens pas la même chose avec...

En fait la différence ne vient pas de lm et glm, mais de anova et summary. On peut vérifier qu'on obtient la même chose avec anova(lm()) et anova(glm()) d'une part, et avec summary(lm()) et summary(glm()) d'autre part.
Il semble qu'avec summary l'ordre n'a plus d'importance parce que chaque effet est testé conditionnellement aux autres (mais il faudrait se plonger dans la doc et vérifier, je n'ai pas le temps), ce qui correspondrait aux sommes des carrés de type III de SAS.

[Logez Maxime]

Bonjour,

La différence que Nico observe provient de control1.2 et de la modalité qui est prise en compte dans le modèle. Normalement il devrait prendre en compte pour l'interaction le fait d'être à la fois 2 de control1 et 2 de control2 donc la modalité 2.2 or dans le modèle de Nico ce n'est pas cette modalité qu'il prend mais la modalité 1.2 :

Code : Tout sélectionner

glm(x~control1+control2+control1.2+etud,data=df,x=TRUE)$x
   (Intercept) control12 control22 control1.21.2 control1.22.1 control1.22.2 etud2 etud3
1            1         0         0             0             0             0     0     0
2            1         0         0             0             0             0     1     0
3            1         0         0             0             0             0     0     1
4            1         0         0             0             0             0     0     1
5            1         0         1             1             0             0     0     0
6            1         0         1             1             0             0     1     0
7            1         1         0             0             1             0     1     0
8            1         1         0             0             1             0     1     0
9            1         1         0             0             1             0     0     1
10           1         1         1             0             0             1     0     0
11           1         1         1             0             0             1     0     1
12           1         1         1             0             0             1     0     1


Iil ne se sert pas de la bonne indicatrice de l'interaction car si on regarde le modèle suivant :

Code : Tout sélectionner

glm(terms(x ~ control1 * control2 + etud, keep.order=T), data=df,x=TRUE)$x
   (Intercept) control12 control22 control12:control22 etud2 etud3
1            1         0         0                   0     0     0
2            1         0         0                   0     1     0
3            1         0         0                   0     0     1
4            1         0         0                   0     0     1
5            1         0         1                   0     0     0
6            1         0         1                   0     1     0
7            1         1         0                   0     1     0
8            1         1         0                   0     1     0
9            1         1         0                   0     0     1
10           1         1         1                   1     0     0
11           1         1         1                   1     0     1
12           1         1         1                   1     0     1


L'indicatrice de l'interaction est control12:control22, soit la modalité 2.2, alors que dans le premier modèle il va se servir de control1.21.2 soit la modalité 1.2.
On peut revenir au même résultat en modifiant les contrasts de l'interaction :

Code : Tout sélectionner

with(df,control1.3 <- control1.2)
contrasts(df$control1.3) <- matrix(c(rep(0,11),1),4)
contrasts(df$control1.3)
    [,1] [,2] [,3]
1.1    0    0    0
1.2    0    0    0
2.1    0    0    0
2.2    0    0    1
glm(x~control1+control2+control1.3+etud,data=df,x=TRUE)$x
   (Intercept) control12 control22 control1.31 control1.32 control1.33 etud2 etud3
1            1         0         0           0           0           0     0     0
2            1         0         0           0           0           0     1     0
3            1         0         0           0           0           0     0     1
4            1         0         0           0           0           0     0     1
5            1         0         1           0           0           0     0     0
6            1         0         1           0           0           0     1     0
7            1         1         0           0           0           0     1     0
8            1         1         0           0           0           0     1     0
9            1         1         0           0           0           0     0     1
10           1         1         1           0           0           1     0     0
11           1         1         1           0           0           1     0     1
12           1         1         1           0           0           1     0     1

# les différents summary :
summary(glm(terms(x ~ control1 * control2 + etud, keep.order=T), data=df))

Call:
glm(formula = terms(x ~ control1 * control2 + etud, keep.order = T),
    data = df)

Deviance Residuals:
     Min        1Q    Median        3Q       Max
-1.04107  -0.19848   0.06375   0.28777   0.54893

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)
(Intercept)          10.1404     0.4349  23.318 4.08e-07 ***
control12            -6.1743     0.4854 -12.720 1.45e-05 ***
control22            -6.5829     0.5455 -12.067 1.97e-05 ***
control12:control22   6.6189     0.8126   8.145 0.000184 ***
etud2                 0.8750     0.5103   1.715 0.137234
etud3                 0.8518     0.4625   1.842 0.115110
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 0.3472179)

    Null deviance: 108.8580  on 11  degrees of freedom
Residual deviance:   2.0833  on  6  degrees of freedom
AIC: 27.043

Number of Fisher Scoring iterations: 2

summary(glm(x~control1+control2+control1.2+etud,data=df))

Call:
glm(formula = x ~ control1 + control2 + control1.2 + etud, data = df)

Deviance Residuals:
     Min        1Q    Median        3Q       Max
-1.04107  -0.19848   0.06375   0.28777   0.54893

Coefficients: (2 not defined because of singularities)
              Estimate Std. Error t value Pr(>|t|)
(Intercept)   10.14036    0.43487  23.318 4.08e-07 ***
control12     -6.17429    0.48540 -12.720 1.45e-05 ***
control22      0.03607    0.55957   0.064 0.950696
control1.21.2 -6.61893    0.81261  -8.145 0.000184 ***
control1.22.1       NA         NA      NA       NA
control1.22.2       NA         NA      NA       NA
etud2          0.87500    0.51031   1.715 0.137234
etud3          0.85179    0.46251   1.842 0.115110
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 0.3472179)

    Null deviance: 108.8580  on 11  degrees of freedom
Residual deviance:   2.0833  on  6  degrees of freedom
AIC: 27.043

Number of Fisher Scoring iterations: 2

summary(glm(x~control1+control2+control1.3+etud,data=df))

Call:
glm(formula = x ~ control1 + control2 + control1.3 + etud, data = df)

Deviance Residuals:
     Min        1Q    Median        3Q       Max
-1.04107  -0.19848   0.06375   0.28777   0.54893

Coefficients: (2 not defined because of singularities)
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  10.1404     0.4349  23.318 4.08e-07 ***
control12    -6.1743     0.4854 -12.720 1.45e-05 ***
control22    -6.5829     0.5455 -12.067 1.97e-05 ***
control1.31       NA         NA      NA       NA
control1.32       NA         NA      NA       NA
control1.33   6.6189     0.8126   8.145 0.000184 ***
etud2         0.8750     0.5103   1.715 0.137234
etud3         0.8518     0.4625   1.842 0.115110
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 0.3472179)

    Null deviance: 108.8580  on 11  degrees of freedom
Residual deviance:   2.0833  on  6  degrees of freedom
AIC: 27.043

Number of Fisher Scoring iterations: 2


Le dernier modèle avec modification de contrasts étant bien égal au premier avec keep.order.

Maxime

[Nicolas Péru]

Donc lorsque l'on veut bricoler une interaction, il vaut mieux jeter un oeil aux contrastes...

[Matthieu Lesnoff]

Si keep.order n'avait pas existé, il aurait fallu l'inventer.

ou dans le cas particulier de ton exemple utiliser drop 1 (qui teste chaque composante déclarée ds le modèle après l'ajustement des autres mais en respectant la 'marginalité') :

Code : Tout sélectionner


> fm1 <- lm(formula = x ~ control1 * control2 + etud, data = df)
> fm2 <- lm(formula = terms(x = x ~ control1 * control2 + etud, keep.order = T), data = df)

> drop1(fm1, test = "F")
Single term deletions

Model:
x ~ control1 * control2 + etud
                  Df Sum of Sq     RSS     AIC F value     Pr(F)   
<none>                          2.0833 -9.0114                     
etud               2    1.4213  3.5046 -6.7699  2.0467 0.2100576   
control1:control2  1   23.0364 25.1197 18.8649 66.3455 0.0001841 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> anova(fm2)
Analysis of Variance Table

Response: x
                  Df Sum Sq Mean Sq  F value    Pr(>F)   
control1           1 43.777  43.777 126.0799 2.981e-05 ***
control2           1 31.507  31.507  90.7425 7.635e-05 ***
control1:control2  1 30.069  30.069  86.5991 8.716e-05 ***
etud               2  1.421   0.711   2.0467    0.2101   
Residuals          6  2.083   0.347                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



Idem pour glm :

Code : Tout sélectionner


> fm3 <- glm(formula = x ~ control1 * control2 + etud, data = df, family = gaussian)
> fm4 <- glm(formula = terms(x = x ~ control1 * control2 + etud, keep.order = T), data = df, family = gaussian)

> drop1(fm3, test = "F")
Single term deletions

Model:
x ~ control1 * control2 + etud
                  Df Deviance    AIC F value     Pr(F)   
<none>                  2.083 27.043                     
etud               2    3.505 29.285  2.0467 0.2100576   
control1:control2  1   25.120 54.919 66.3455 0.0001841 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> anova(fm4, test = "F")
Analysis of Deviance Table

Model: gaussian, link: identity

Response: x

Terms added sequentially (first to last)


                  Df Deviance Resid. Df Resid. Dev        F    Pr(>F)   
NULL                                 11    108.858                       
control1           1   43.777        10     65.081 126.0799 2.981e-05 ***
control2           1   31.507         9     33.573  90.7425 7.635e-05 ***
control1:control2  1   30.069         8      3.505  86.5991 8.716e-05 ***
etud               2    1.421         6      2.083   2.0467    0.2101   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



Ci-dessus, le test de 'etud' est un test de 'type III' sous SAS.

Avec summary, chaque *coefficient* du modèle (= chaque composante du vecteur 'b' dans 'y = X b') est testé à zéro après l'ajustement de tous les autres, ceci en effet indépendemment de l'ordre d'entrée des facteurs ou covariables.

Matthieu

[Alejandro Collantes]

Bonjour
Il y a une question similaire a le FAQ de R....

http://cran.r-project.org/doc/FAQ/R-FAQ.html

7.18 Why does the output from anova() depend on the order of factors in the model?
In a model such as ~A+B+A:B, R will report the difference in sums of squares between the models ~1, ~A, ~A+B and ~A+B+A:B. If the model were ~B+A+A:B, R would report differences between ~1, ~B, ~A+B, and ~A+B+A:B . In the first case the sum of squares for A is comparing ~1 and ~A, in the second case it is comparing ~B and ~B+A. In a non-orthogonal design (i.e., most unbalanced designs) these comparisons are (conceptually and numerically) different.

Some packages report instead the sums of squares based on comparing the full model to the models with each factor removed one at a time (the famous `Type III sums of squares' from SAS, for example). These do not depend on the order of factors in the model. The question of which set of sums of squares is the Right Thing provokes low-level holy wars on R-help from time to time.

There is no need to be agitated about the particular sums of squares that R reports. You can compute your favorite sums of squares quite easily. Any two models can be compared with anova(model1, model2), and drop1(model1) will show the sums of squares resulting from dropping single terms.

Cordielment,