Cases omises de summary.glm()

[Nicolas Mazziotta]

Bonjour,

Je poursuis ma découverte des méthodes loglinéaires. J'aimerais à présent évaluer la significativité des effets à l'aide de tests de Wald. Comment faire quand les effets sont considérés comme redondants et n'apparaissent pas dans summary.glm(). Par exemple, soit la table

Code : Tout sélectionner

   MI MM PPD Freq
1   a  w   0   81
2   b  w   0  159
3   a  x   0  251
4   b  x   0  283
5   a  y   0    8
6   b  y   0   20
7   a  z   0   12
8   b  z   0   20
9   a  w   1   13
10  b  w   1   74
11  a  x   1   24
12  b  x   1   77
13  a  y   1    4
14  b  y   1    8
15  a  z   1    2
16  b  z   1    9


Je sélectionne le modèle

Code : Tout sélectionner

Call:
glm(formula = Freq ~ MM + MI + PPD + MM:MI + MM:PPD + MI:PPD,
    family = poisson, data = data.frame(tmp))

Deviance Residuals:
        1          2          3          4          5          6          7
 0.090599  -0.064338   0.070038  -0.065776  -0.637618   0.449021   0.006014
        8          9         10         11         12         13         14
-0.004654  -0.220857   0.094817  -0.222774   0.126876   1.196424  -0.637618
       15         16
-0.014673   0.006946

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  4.38437    0.10714  40.922  < 2e-16 ***
MMx          1.13666    0.12174   9.337  < 2e-16 ***
MMy         -2.08765    0.31298  -6.670 2.55e-11 ***
MMz         -1.90120    0.29162  -6.519 7.06e-11 ***
MIb          0.68964    0.12842   5.370 7.86e-08 ***
PPD1        -1.75878    0.19733  -8.913  < 2e-16 ***
MMx:MIb     -0.56131    0.14856  -3.778 0.000158 ***
MMy:MIb     -0.09273    0.37231  -0.249 0.803308
MMz:MIb     -0.17603    0.35318  -0.498 0.618188
MMx:PPD1    -0.53906    0.16876  -3.194 0.001402 **
MMy:PPD1     0.18407    0.37329   0.493 0.621933
MMz:PPD1    -0.02088    0.37722  -0.055 0.955847
MIb:PPD1     0.97780    0.18646   5.244 1.57e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 1502.8536  on 15  degrees of freedom
Residual deviance:    2.5916  on  3  degrees of freedom
AIC: 110.18

Number of Fisher Scoring iterations: 4


Je peux facilement déduire l'effet et la probabilité de PPD0 à partir de ceux de PPD1, puisque la variable est dichotomique. Mais comment faire concrètement quand la variable n'est pas dichotomique. Par exemple dans le cas de MMw ou de MMw:MIb? J'ai fouillé dans les objets R glm et summary.glm(), mais les effets redondants sont éliminés d'office. Comment les restituer?

Merci pour toute aide.

[Renaud Lancelot]

Il n'y a pas d'"effet redondant"... Les variables qualitatives - ou facteurs dans R - sont transformées en indicatrices avant les calculs. Cette transformation est effectuée selon des règles définies par les contrastes. Par défaut, dans R, le contraste utilisé pour les facteurs est "treatment", qui aboutit à prendre comme modalité de référence le premier niveau du facteur. En fin de compte, les variables à expliquer et explicatives sont transformées en une matrice (avec la fct model.matrix) utilisée dans l'analyse:

Code : Tout sélectionner

> x <- factor(c("a", "b"))
> x
[1] a b
Levels: a b
> options()$contrasts
        unordered           ordered
"contr.treatment"      "contr.poly"
> x <- factor(c("a", "b"))
> x
[1] a b
Levels: a b
> options()$contrasts
        unordered           ordered
"contr.treatment"      "contr.poly"
> contrasts(x)
  b
a 0
b 1
> model.matrix(~ x)
  (Intercept) xb
1           1  0
2           1  1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$x
[1] "contr.treatment"


Le niveau de base du facteur se retrouve ainsi dans la constante du modèle. Quand il y a plusieurs facteurs, la constante s'interprète comme la moyenne de pop correspondant aux niveaux de base des différents facteurs:

Code : Tout sélectionner

> z <- factor(c("c", "d"))
> model.matrix(~ x + z)
  (Intercept) xb zd
1           1  0  0
2           1  1  1
attr(,"assign")
[1] 0 1 2
attr(,"contrasts")
attr(,"contrasts")$x
[1] "contr.treatment"

attr(,"contrasts")$z
[1] "contr.treatment"

Dans un modèle comportant des facteurs à plusieurs modalités (> 2),
* vous pouvez faire des tests d'égalité d'un ou plusieurs coefs à des valeurs arbitraires (tests de Wald: cf précédents échanges).
* Si vous voulez tester globalement l'effet d'un facteur ou d'une interaction, ajuster les modèles avec et sans ces effet ou interaction et les comparer avec un test du rapport des vraisemblances.

Code : Tout sélectionner

> library(MASS)
> fm1 <- glm(Postwt ~ Prewt + offset(Prewt), family = gaussian, data = anorexia)
> fm2 <- glm(Postwt ~ Prewt + Treat + offset(Prewt), family = gaussian, data = anorexia)
> # test effet Treat
> anova(fm1, fm2, test = "Chisq")
Analysis of Deviance Table

Model 1: Postwt ~ Prewt + offset(Prewt)
Model 2: Postwt ~ Prewt + Treat + offset(Prewt)
  Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1        70     4077.5                     
2        68     3311.3  2    766.3 0.0003828


[Nicolas Péru]

Le fait que la première modalité d'un facteur (première lettre de l'alphabet s'il est codé en lettre et premier chiffre dans l'ordre si c'est un nombre) n'est pas dû au hasard ou simplement au bon vouloir des développeurs de R.
En fait les fonction lm et autre glm sont construites pour tester des plans d'expérience où il y a un "tube témoin" comme en chimie par exemple. Par défaut, ce témoin est la première modalité du facteur. Donc la significativité donnée par R de l'effet de la modalité du facteur n'est pas directement celle sur la variable expliquée.
Je traduis avec l'exemple :
pour la variable MM, c'est la modalité w qui sert de témoin. On voit que pour la modalité x, R donne une p-value très faible donc le paramètre est signifiactivement différent de zéro. De fait, l'effet de la modalité x par rapport à la modalité w prise seule est significatif, de fait l'ajout de la modalité x apporte quelque chose au modèle par rapport à prendre uniquement la modalité w.
Pour tester l'effet de la variable catégorielle une simple anova suffit. Si tu veux plus de détail quant aux effets de chaque modalité du facteur la fonction summary est plus adaptée.
Donc tu devras toujours avoir une modalité témoin (à ma connaissance en tout cas). Modifier la matrice des contrastes te permet de changer de témoin et voir l'effet sur ton modèle : est ce qu'il ya une modalité qui apporte plus au modèle que les autres ?
Dans le cas d'étude avec des variables catégorielles, le seul moyen de voir apparaitre la première modalité dans le summary c'set d'enlever l'intercept. Par contre, je ne saurais te répondre quant au bien fondé de cette manipulation (variable~-1+facteur1+facteur2...)


enfin pour ce qu'il s'agit de manipuler les contrastes je te recommande cette fiche très bien faite par Daniel Chessel : http://biom3.univ-lyon1.fr/R/fichestd/tdr334.pdf
Nicolas

[Tillard]

Bonjour

Code : Tout sélectionner

Modifier la matrice des contrastes te permet de changer de témoin et voir l'effet sur ton modèle : est ce qu'il ya une modalité qui apporte plus au modèle que les autres ?


voir aussi la fonction relevel() qui facilite le changement de la modalité de référence d'un facteur; on peut difficilement faire plus simple;
cordialement

[Nicolas Mazziotta]

Si je comprends bien, l'estimate et la SE de la ligne

Code : Tout sélectionner

MMx:MIb     -0.56131    0.14856  -3.778 0.000158 ***


Dans mon exemple permet d'estimer l'oddsratio de la table 2x2 croisant MM (modalités w et x uniquement) et MI (modalités a et b). Je peux donc estimer l'intervalle de confiance comme suit

Code : Tout sélectionner

>1/exp(-.56+1.96*.15)
[1] 1.304735
> 1/exp(-.56-1.96*.15)
[1] 2.349024


et affirmer que l'association de MMx à MMa a entre 1.3 et 2.3 fois plus de "chances" de se réaliser que celle de MMx à MMb.

Par contre, je ne peux rien dire de l'opposition entre MMx et MMy, parce que la référence est MMw. Il faudrait pour ce faire changer l'ordre des niveaux. Je dois donc créer plusieurs objets glm équivalents (même AIC) pour évaluer tous les oddsratios.

[Renaud Lancelot]

Euh, je ne suis pas sûr que l'estimation des odds ratios soit pertinente dans un cas comme celui-ci. Quelle est la question à laquelle vous voulez répondre ?

Renaud

[Nicolas Mazziotta]

Je vais essayer d'être clair.

Les procédures de sélection de modèle m'ont mené à choisir le modèle MI+MM+PPD+MI:PPD+MI:MM+MM:PPD pour ma table tridimensionnelle.

À présent, la question de savoir si l'une ou l'autre de ces variables est plus ou moins impliquée ne m'intéresse plus (la question est résolue). Je veux savoir, pour chaque association de deux variables, dans quel sens vont les associations:

Quand les deux variables sont dichotomiques (comme dans l'échange viewtopic.php?t=290), c'est très facile, mais dans le cas d'un tableau IxJ, je ne vois pas concrètement. Par exemple, j'ai le tableau de deux colonnes et quatre lignes MM:MI. Les résultats des tests de Wald donnés par summary.glm() sont:

Code : Tout sélectionner

MMx:MIb     -0.56131    0.14856  -3.778 0.000158 ***
MMy:MIb     -0.09273    0.37231  -0.249 0.803308
MMz:MIb     -0.17603    0.35318  -0.498 0.618188


Avec ça, je peux juste me faire une idée des attractions comme suit:

Code : Tout sélectionner

    MIa MIb
MMw  ?   ?
MMx  ?   -
MMy  ?   -(non signif.)
MMz  ?   -(non signif.)


Mais j'aimerais décomposer la table de contingence comme je l'ai fait dans l'échange viewtopic.php?t=326.

Il se peut que quelque chose m'ait échappé et que la réponse ait déjà été donnée.

J'espère avoir été clair. Merci encore.

[Nicolas Mazziotta]

Je vais approfondir mes connaissances en la matière avant de poster à nouveau. Merci pour votre aide jq ici.

À bientôt.

[Logez Maxime]

Bonjour,

Je vais encore faire le rabat joie mais ici lorsque tu regardes les tets de wald, selon quel levels du facteur est pris comme référence dans ton modèle tu n'auras pas les mêmes conclusions :

Code : Tout sélectionner

test
   MI MM PPD Freq
1   a  w   0   81
2   b  w   0  159
3   a  x   0  251
4   b  x   0  283
5   a  y   0    8
6   b  y   0   20
7   a  z   0   12
8   b  z   0   20
9   a  w   1   13
10  b  w   1   74
11  a  x   1   24
12  b  x   1   77
13  a  y   1    4
14  b  y   1    8
15  a  z   1    2
16  b  z   1    9

# premier modèle avec les premières modalités par défaut :
glm1 <- glm(formula = Freq ~ MM + MI + PPD + MM:MI + MM:PPD + MI:PPD, family = poisson, data = data.frame(test))

summary(glm1)
Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  4.38437    0.10714  40.922  < 2e-16 ***
MMx          1.13666    0.12174   9.337  < 2e-16 ***
MMy         -2.08765    0.31298  -6.670 2.55e-11 ***
MMz         -1.90120    0.29162  -6.519 7.06e-11 ***
MIb          0.68964    0.12842   5.370 7.86e-08 ***
PPD1        -1.75878    0.19733  -8.913  < 2e-16 ***
MMx:MIb     -0.56131    0.14856  -3.778 0.000158 ***
MMy:MIb     -0.09273    0.37231  -0.249 0.803308   
MMz:MIb     -0.17603    0.35318  -0.498 0.618188   
MMx:PPD1    -0.53906    0.16876  -3.194 0.001402 **
MMy:PPD1     0.18407    0.37329   0.493 0.621933   
MMz:PPD1    -0.02088    0.37722  -0.055 0.955847   
MIb:PPD1     0.97780    0.18646   5.244 1.57e-07 ***

# ici interaction entre MMx et MIb significative

# le même modèle mais en fixant la modalité z de MM comme référence à la place de w
contrasts(test[,2])
  x y z
w 0 0 0
x 1 0 0
y 0 1 0
z 0 0 1
contrasts(test[,2]) <- contr.treatment(levels(test[,2]),4)
 
contrasts(test[,2])
  w x y
w 1 0 0
x 0 1 0
y 0 0 1
z 0 0 0

glm2 <- glm(formula = Freq ~ MM + MI + PPD + MM:MI + MM:PPD + MI:PPD,family = poisson, data = data.frame(test))
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  2.48317    0.27295   9.098  < 2e-16 ***
MMw          1.90120    0.29162   6.519 7.06e-11 ***
MMx          3.03786    0.27905  10.886  < 2e-16 ***
MMy         -0.18646    0.40129  -0.465    0.642   
MIb          0.51360    0.33240   1.545    0.122   
PPD1        -1.77967    0.38391  -4.636 3.56e-06 ***
MMw:MIb      0.17603    0.35318   0.498    0.618   
MMx:MIb     -0.38528    0.34067  -1.131    0.258   
MMy:MIb      0.08330    0.48227   0.173    0.863   
MMw:PPD1     0.02088    0.37722   0.055    0.956   
MMx:PPD1    -0.51818    0.37194  -1.393    0.164   
MMy:PPD1     0.20496    0.49933   0.410    0.681   
MIb:PPD1     0.97780    0.18646   5.244 1.57e-07 ***

# ici interaction entre MMx et MIb non significative

# par contre :
anova(glm1,test="Chisq")
Analysis of Deviance Table

Model: poisson, link: log

Response: Freq

Terms added sequentially (first to last)


       Df Deviance Resid. Df Resid. Dev  P(>|Chi|)
NULL                      15    1502.85           
MM      3   969.46        12     533.39 7.577e-210
MI      1    62.86        11     470.53  2.222e-15
PPD     1   397.32        10      73.22  2.115e-88
MM:MI   3    21.49         7      51.73  8.331e-05
MM:PPD  3    18.34         4      33.38  3.734e-04
MI:PPD  1    30.79         3       2.59  2.871e-08

anova(glm2,test="Chisq")
Analysis of Deviance Table

Model: poisson, link: log

Response: Freq

Terms added sequentially (first to last)


       Df Deviance Resid. Df Resid. Dev  P(>|Chi|)
NULL                      15    1502.85           
MM      3   969.46        12     533.39 7.577e-210
MI      1    62.86        11     470.53  2.222e-15
PPD     1   397.32        10      73.22  2.115e-88
MM:MI   3    21.49         7      51.73  8.331e-05
MM:PPD  3    18.34         4      33.38  3.734e-04
MI:PPD  1    30.79         3       2.59  2.871e-08


Conclusions : le test de wald est intéressant mais dépend des modalités prises comme référence, et ne permet pas ici de conclure si une modalité d'un facteur est liée à une autre modalité d'un deuxième facteur. La seule info je pense que tu peux tirer de ces modèles est une information générale qui te dit si les interactions sont significatives ou non (via les anova).

Maxime

[Nicolas Mazziotta]

Dommage.
Merci tout de même