[Réglé] Modèles log-linéaires: ai-je compris?

Postez ici vos questions, réponses, commentaires ou suggestions - Les sujets seront ultérieurement répartis dans les archives par les modérateurs

Modérateur : Groupe des modérateurs

Nicolas Mazziotta
Messages : 18
Enregistré le : 02 Jan 2007, 09:43

[Réglé] Modèles log-linéaires: ai-je compris?

Messagepar Nicolas Mazziotta » 17 Jan 2007, 07:43

Bonjour à tous,

Tout d'abord, veuillez m'excuser si je profère des insanités: j'ai une formation plutôt littéraire et je me suis mis aux stats il y a peu et j'essaye de me servir de R.

Je voulais me servir d'une modélisation log-linéaire pour étudier la table tridimensionnelle suivante:

Code : Tout sélectionner


      E     P     H  Fq
1 FALSE  TRUE  TRUE  52
2  TRUE  TRUE  TRUE   2
3 FALSE FALSE  TRUE 802
4  TRUE FALSE  TRUE  26
5 FALSE  TRUE FALSE  19
6  TRUE  TRUE FALSE  39
7 FALSE FALSE FALSE 195
8  TRUE FALSE FALSE 196




Je construit les modèles et je les teste:

Code : Tout sélectionner

> m1<-glm(Fq ~ E*H*P ,poisson,tt)
> anova(m1, test="Chisq")
Analysis of Deviance Table

Model: poisson, link: log

Response: Fq

Terms added sequentially (first to last)


      Df Deviance Resid. Df Resid. Dev  P(>|Chi|)
NULL                      7    2218.26
E      1   522.01         6    1696.26 1.548e-115
H      1   143.46         5    1552.80  4.664e-33
P      1  1076.42         4     476.38 4.412e-236
E:H    1   453.39         3      22.99 1.322e-100
E:P    1    18.98         2       4.02  1.324e-05
H:P    1     3.51         1       0.51       0.06
E:H:P  1     0.51         0  2.776e-14       0.48


Si je comprends bien, avec un seuil alpha traditionnel, je dois rejeter le modèle saturé et considérer que le meilleur modèle est: Freq ~ H:E+P:E

En d'autres termes, les interactions sont groupées deux à deux: chaque variable influence une des deux autres à la fois. Mais il n'y a pas d'interaction significative entre H et P.

Logez Maxime
Messages : 3138
Enregistré le : 26 Sep 2006, 11:35

Messagepar Logez Maxime » 17 Jan 2007, 15:39

Bonjour,


Par ton modèle que cherches-tu à mettre en évidence ?

Maxime

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 18 Jan 2007, 07:36

Deux remarques:

1. Il est dangereux d'avoir des noms de variables avec des points dans des formules. Le point "." a en effet une signification particulière dans les formules: voir par exemple "An introduction to R" (livré avec R), p.54 pour la version de R 2.4.1

2. Il est préférable de mettre les données sous formes d'un data.frame, ce qui a le mérite de permettre une représentation explicite des données par rapport à la formule. Dans votre exemple, j'ai ainsi du mal à voir ce qui est modélisé.

Renaud

Nicolas Mazziotta
Messages : 18
Enregistré le : 02 Jan 2007, 09:43

Messagepar Nicolas Mazziotta » 18 Jan 2007, 17:26

J'ai édité le post pour qu'il soit plus lisible.

Je travaille sur la ponctuation médiévale. Je cherche à déterminer comment interagissent la ponctuation initiale des phrases (P), la présence d'une majuscule (H) et la présence du mot _et_ au début de la phrase (E).

Si j'ai bien compris, le modèle me montre que le mot _et_ interagit avec la ponctuation et avec la majuscule, mais que ces deux derniers ne s'influencent pas mutuellement

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 18 Jan 2007, 18:55

Il vaut mieux utiliser la sortie de la fonction summary (tests de Wald) ou des tests du rapport des vraisemblances (LRT) pour juger de l'effet d'un facteur, la seconde option étant préférable pour des raisons statistiques (pb de l'effet Hauck-Donner avec le test de Wald).

Code : Tout sélectionner

> m1 <- glm(Fq ~ E * H * P , family= poisson, data = tt)
> summary(m1)

Call:
glm(formula = Fq ~ E * H * P, family = poisson, data = tt)

Deviance Residuals:
[1]  0  0  0  0  0  0  0  0

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)   
(Intercept)        5.273000   0.071611  73.633   <2e-16 ***
ETRUE              0.005115   0.101145   0.051   0.9597   
HTRUE              1.414109   0.079844  17.711   <2e-16 ***
PTRUE             -2.328561   0.240333  -9.689   <2e-16 ***
ETRUE:HTRUE       -3.434127   0.223470 -15.367   <2e-16 ***
ETRUE:PTRUE        0.714008   0.297494   2.400   0.0164 * 
HTRUE:PTRUE       -0.407304   0.279710  -1.456   0.1453   
ETRUE:HTRUE:PTRUE -0.543092   0.804638  -0.675   0.4997   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance:  2.2183e+03  on 7  degrees of freedom
Residual deviance: -6.4393e-14  on 0  degrees of freedom
AIC: 62.559

Number of Fisher Scoring iterations: 3


Par exemple, le LRT pour l'interaction d'ordre 2 (3 voies):

Code : Tout sélectionner

> m2 <- update(m1, formula = . ~ . - E:H:P)
>
> anova(m1, m2, test = "Chisq")
Analysis of Deviance Table

Model 1: Fq ~ E * H * P
Model 2: Fq ~ E + H + P + E:H + E:P + H:P
  Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1         0 -6.439e-14                     
2         1    0.50620 -1 -0.50620   0.47679


Dans ce cas, on a une bonne correspondance entre le test de Wald et le LRT.

Ensuite, il est préférable d'utiliser une procédure appropriée pour sélectionner le "meilleur" modèle. Les comparaisons multiples de modèles pris 2 à 2 sont à proscrire (pas de théorie statistique correcte pour les comparaisons). La fonction "step" offre une procédure pas-à-pas basée sur le critère d'information d'Akaike (AIC) ou d'autres critères: voir l'aide pour des explications.

Code : Tout sélectionner

> m3 <- step(m1, scope = list(lower = . ~ 1, upper = . ~ .), dir = "both")
Start:  AIC= 62.56
 Fq ~ E * H * P

        Df   Deviance    AIC
- E:H:P  1      0.506 61.065
<none>     -6.439e-14 62.559

Step:  AIC= 61.07
 Fq ~ E + H + P + E:H + E:P + H:P

        Df   Deviance    AIC
<none>           0.51  61.07
+ E:H:P  1 -6.439e-14  62.56
- H:P    1       4.02  62.58
- E:P    1       6.11  64.67
- E:H    1     440.52 499.08
> summary(m3)

Call:
glm(formula = Fq ~ E + H + P + E:H + E:P + H:P, family = poisson,
    data = tt)

Deviance Residuals:
       1         2         3         4         5         6         7         8 
 0.13673  -0.60393  -0.03458   0.19461  -0.22102   0.15822   0.07028  -0.06987 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  5.26796    0.07141  73.773   <2e-16 ***
ETRUE        0.01514    0.10007   0.151   0.8798   
HTRUE        1.42037    0.07946  17.874   <2e-16 ***
PTRUE       -2.27324    0.22112 -10.280   <2e-16 ***
ETRUE:HTRUE -3.48378    0.21439 -16.249   <2e-16 ***
ETRUE:PTRUE  0.62826    0.26657   2.357   0.0184 * 
HTRUE:PTRUE -0.48286    0.25313  -1.908   0.0564 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 2218.2642  on 7  degrees of freedom
Residual deviance:    0.5062  on 1  degrees of freedom
AIC: 61.065

Number of Fisher Scoring iterations: 4


Noter que les résultats peuvent être très différents de ceux obtenus avec anova, de la manière dont vous l'avez utilisée: une anova séquentielle est réalisée, les effets étant ajoutés les uns après les autres.

Pour l'interprétation, on peut dire que la ponctuation initiale dépend de la présence d'une majuscule et de la présence d'un "et" en début de phrase (et réciproquement), mais que la ponctuation initiale associée à la présence d'une majuscule ne dépend pas de la présence de "et" en début de phrase (et de même avec les autres combinaisons: P et E ne dépend pas de H, E et H ne dépend pas de P).

Renaud

Nicolas Mazziotta
Messages : 18
Enregistré le : 02 Jan 2007, 09:43

Messagepar Nicolas Mazziotta » 18 Jan 2007, 20:12

Merci pour cette réponse rapide et riche.

Nicolas Mazziotta
Messages : 18
Enregistré le : 02 Jan 2007, 09:43

Coefficients et comparaisons multiples

Messagepar Nicolas Mazziotta » 16 Fév 2007, 09:27

Je me permets de revenir sur la question.
Le modèle choisi permet de de faire une idée de l'importance des différentes valeurs des variables à l'aide des probabilités des coefficients.

Ma question est la suivante: chacun des tests évaluant les probabilité est indépendant, mais le fait de les considérer tous ensemble mène à un problème de comparaison multiple. Dans ce cas, pour évaluer en toute rigueur l'intérêt de chaque variable, il faudrait appliquer une correction de Bonferroni ou Sidak (pour tests indépendants). Suis-je dans le bon?

S'il faut effectivement opérer un ajustement, la question suivante est: combien de tests considérer? Faut-il uniquement considérer les tests non redondants repris par summary() ou tous les tests possibles?

Merci pour toute aide/tout commentaire.

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 16 Fév 2007, 18:13

La question de la comparaison de modèles (ou des tests multiples, de manière équivalente) est très vaste et débattue. Pour ma part, j'ai recours aux critères d'information pour lesquels une très brève introduction est disponible sur ce site.

En bref, il s'agit d'établir la liste des modèles plausibles, de les estimer et de les comparer avec un critère d'information (AIC, AICc,...).

Voir les références citées dans la fiche.

Renaud

Nicolas Mazziotta
Messages : 18
Enregistré le : 02 Jan 2007, 09:43

Messagepar Nicolas Mazziotta » 17 Fév 2007, 11:08

Merci pour ce lien.
J'avais bien compris que l'AIC (ou le BIC) permettait de "comparer" les modèles.

Ceci dit, je me posais la question des valeurs de probabilité des coefficients. Dans la sortie de summary les coefficients (ici, ETRUE,PTRUE, etc.) sont associés à une probabilité liée à z: est-elle utilisable en l'état ou le recours à un ajustement est-il préférable.

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 17 Fév 2007, 17:35

Pour un modèle donné, vous pouvez utiliser ces valeurs de P qui correspondent à des tests de Wald, avec leurs avantages et leurs inconvénients (problème de l'effet Hauck-Donner). Les méthodes citées dans la fiche (Bayesian model averaging, inférence multi-modèle,...) permettent de tenir compte de l'incertitude sur la sélection des modèles pour juger de l'importance de telle ou telle variable.

Renaud

Nicolas Mazziotta
Messages : 18
Enregistré le : 02 Jan 2007, 09:43

Messagepar Nicolas Mazziotta » 18 Fév 2007, 15:11

Merci.

Finalement, quand toutes les variables sont dichotomiques, j'emploie la fonction dropterm(MASS) pour obtenir les LRT des associations. La sortie synthétise les résultats des anova du modèle choisi comparé aux modèles amputés d'un terme (ici, trois comparaisons)

Code : Tout sélectionner

>dropterm(m3, test="Chisq")
Single term deletions

Model:
Fq ~ E + H + P + E:H + E:P + H:P
       Df Deviance    AIC    LRT Pr(Chi)
<none>        0.51  61.07
E:H     1   440.52 499.08 440.01 < 2e-16 ***
E:P     1     6.11  64.67   5.60 0.01792 *
H:P     1     4.02  62.58   3.51 0.06089 .


Malheureusement (désolé d'insister), ce sont trois anova différentes qui sont effectuées. Les probabilités associées aux LRT sont donc multiples. Si je veux rester orthodoxe, je dois donc multiplier par trois (Bonferroni) la colonne Pr(Chi)? Sinon, comment se fait-il que la procédure convienne?

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 18 Fév 2007, 22:23

Trois comparaisons, c'est encore raisonnable... Au-delà, employer une technique de sélection des modèles reposant sur l'aic ou autres critères, et tenir compte de l'incertitude sur la sélection des modèles pour votre inférence (valeur de P). Voir les références que je vous ai indiquées.

Renaud

Nicolas Mazziotta
Messages : 18
Enregistré le : 02 Jan 2007, 09:43

Messagepar Nicolas Mazziotta » 19 Fév 2007, 08:33

D'accord. Le package BMA m'a l'air de convenir. J'ai d'autres questions concernant les odds ratios (excusez ma naïveté) pour lesquelles je vais ouvrir un nouveau topic.

Merci pour votre aide précieuse.


Retourner vers « Questions en cours »

Qui est en ligne

Utilisateurs parcourant ce forum : Google [Bot] et 1 invité